Vectorización de Matlab: índices de fila de matriz no cero a la celda

10

Estoy trabajando con Matlab.

Tengo una matriz cuadrada binaria. Para cada fila, hay una o más entradas de 1. Quiero pasar por cada fila de esta matriz y devolver el índice de esos 1s y almacenarlos en la entrada de una celda.

Me preguntaba si hay una manera de hacer esto sin recorrer todas las filas de esta matriz, ya que el ciclo es realmente lento en Matlab.

Por ejemplo, mi matriz

M = 0 1 0
    1 0 1
    1 1 1 

Luego, eventualmente, quiero algo como

A = [2]
    [1,3]
    [1,2,3]

Entonces Aes una célula.

¿Hay alguna manera de lograr este objetivo sin usar for loop, con el objetivo de calcular el resultado más rápidamente?

ftxx
fuente
¿Desea que el resultado sea rápido o desea que el resultado evite forbucles? Para este problema, con las versiones modernas de MATLAB, sospecho que un forbucle es la solución más rápida. Si tiene un problema de rendimiento, sospecho que está buscando la solución en el lugar equivocado en base a consejos desactualizados.
Será
@Quiero que los resultados sean rápidos. Mi matriz es muy grande. El tiempo de ejecución es de aproximadamente 30 segundos en mi computadora usando for loop. Quiero saber si hay algunas operaciones de vectorización inteligentes o, mapReduce, etc. que pueden aumentar la velocidad.
ftxx
1
Sospecho que no puedes. La vectorización funciona en vectores y matrices descritos con precisión, pero su resultado permite vectores de diferentes longitudes. Por lo tanto, mi suposición es que siempre tendrá algún bucle explícito o algún bucle disfrazado cellfun.
HansHirse
@ftxx ¿qué tan grande? ¿Y cuántos 1s en una fila típica? No esperaría que un findciclo tome algo cercano a 30 años para algo lo suficientemente pequeño como para caber en la memoria física.
Será
@ftxx Por favor, vea mi respuesta actualizada, la he editado desde que fue aceptada con una mejora de rendimiento menor
Wolfie

Respuestas:

11

Al final de esta respuesta hay un código de evaluación comparativa, ya que usted aclaró que está interesado en el rendimiento en lugar de evitar arbitrariamente los forbucles.

De hecho, creo que los forbucles son probablemente la opción más eficaz aquí. Desde que se introdujo el "nuevo" (2015b) motor JIT ( fuente ), los forbucles no son inherentemente lentos, de hecho, están optimizados internamente.

Puede ver desde el punto de referencia que la mat2cellopción ofrecida por ThomasIsCoding aquí es muy lenta ...

Comparación 1

Si nos deshacemos de esa línea para splitapplyaclarar la escala, entonces mi método es bastante lento, la opción accumarray de obchardon es un poco mejor, pero las opciones más rápidas (y comparables) están usando arrayfun(como también lo sugiere Thomas) o un forbucle. Tenga en cuenta que arrayfunes básicamente un forbucle disfrazado para la mayoría de los casos de uso, ¡así que no es un lazo sorprendente!

Comparación 2

Le recomendaría que use un forbucle para aumentar la legibilidad del código y el mejor rendimiento.

Editar :

Si suponemos que el bucle es el enfoque más rápido, podemos hacer algunas optimizaciones en torno al findcomando.

Específicamente

  • Hacer Mlógico Como se muestra en el gráfico a continuación, esto puede ser más rápido para relativamente pequeño M, pero más lento con el intercambio de conversión de tipo para grande M.

  • Use un lógico Mpara indexar una matriz en 1:size(M,2)lugar de usar find. Esto evita la parte más lenta del bucle (el findcomando) y supera la sobrecarga de conversión de tipo, por lo que es la opción más rápida.

Aquí está mi recomendación para el mejor rendimiento:

function A = f_forlooplogicalindexing( M )
    M = logical(M);
    k = 1:size(M,2);
    N = size(M,1);
    A = cell(N,1);
    for r = 1:N
        A{r} = k(M(r,:));
    end
end

He agregado esto al punto de referencia a continuación, aquí está la comparación de los enfoques de estilo de bucle:

Comparación 3

Código de evaluación comparativa:

rng(904); % Gives OP example for randi([0,1],3)
p = 2:12; 
T = NaN( numel(p), 7 );
for ii = p
    N = 2^ii;
    M = randi([0,1],N);

    fprintf( 'N = 2^%.0f = %.0f\n', log2(N), N );

    f1 = @()f_arrayfun( M );
    f2 = @()f_mat2cell( M );
    f3 = @()f_accumarray( M );
    f4 = @()f_splitapply( M );
    f5 = @()f_forloop( M );
    f6 = @()f_forlooplogical( M );
    f7 = @()f_forlooplogicalindexing( M );

    T(ii, 1) = timeit( f1 ); 
    T(ii, 2) = timeit( f2 ); 
    T(ii, 3) = timeit( f3 ); 
    T(ii, 4) = timeit( f4 );  
    T(ii, 5) = timeit( f5 );
    T(ii, 6) = timeit( f6 );
    T(ii, 7) = timeit( f7 );
end

plot( (2.^p).', T(2:end,:) );
legend( {'arrayfun','mat2cell','accumarray','splitapply','for loop',...
         'for loop logical', 'for loop logical + indexing'} );
grid on;
xlabel( 'N, where M = random N*N matrix of 1 or 0' );
ylabel( 'Execution time (s)' );

disp( 'Done' );

function A = f_arrayfun( M )
    A = arrayfun(@(r) find(M(r,:)),1:size(M,1),'UniformOutput',false);
end
function A = f_mat2cell( M )
    [i,j] = find(M.');
    A = mat2cell(i,arrayfun(@(r) sum(j==r),min(j):max(j)));
end
function A = f_accumarray( M )
    [val,ind] = ind2sub(size(M),find(M.'));
    A = accumarray(ind,val,[],@(x) {x});
end
function A = f_splitapply( M )
    [r,c] = find(M);
    A = splitapply( @(x) {x}, c, r );
end
function A = f_forloop( M )
    N = size(M,1);
    A = cell(N,1);
    for r = 1:N
        A{r} = find(M(r,:));
    end
end
function A = f_forlooplogical( M )
    M = logical(M);
    N = size(M,1);
    A = cell(N,1);
    for r = 1:N
        A{r} = find(M(r,:));
    end
end
function A = f_forlooplogicalindexing( M )
    M = logical(M);
    k = 1:size(M,2);
    N = size(M,1);
    A = cell(N,1);
    for r = 1:N
        A{r} = k(M(r,:));
    end
end
Wolfie
fuente
1
Ya vi y voté. :-) Todavía estoy esperando a Luis; seguro que tiene algo de magia negra MATLAB para eso.
HansHirse
@Hans Haha, sí, aunque su bolsa habitual de trucos (expansión implícita, indexación inteligente, ...) generalmente mantiene las cosas como matrices, el cuello de botella aquí se resume en celdas
Wolfie
1
Tenga en cuenta que estos tiempos dependen en gran medida de la escasez de M. Si, por ejemplo, solo se llena M = randi([0,20],N) == 20;el 5% de los elementos, entonces el forciclo es, con mucho, el más lento y su arrayfunmétodo gana.
Será
@HansHirse :-) Mi enfoque hubiera sido accumarraysin ind2sub, pero es más lento que el forciclo
Luis Mendo
2

Puede intentar arrayfuncomo a continuación, que barre filas deM

A = arrayfun(@(r) find(M(r,:)),1:size(M,1),'UniformOutput',false)

A =
{
  [1,1] =  2
  [1,2] =

     1   3

  [1,3] =

     1   2   3

}

o (un enfoque más lento por mat2cell)

[i,j] = find(M.');
A = mat2cell(i,arrayfun(@(r) sum(j==r),min(j):max(j)))

A =
{
  [1,1] =  2
  [2,1] =

     1
     3

  [3,1] =

     1
     2
     3

}
ThomasIsCoding
fuente
1
Aunque arrayfunes básicamente un disfraz de bucle, por lo que esto puede fallar en ambos frentes de 1) evitar bucles y 2) ser rápido, como esperaba el OP
Wolfie
2

Editar : agregué un punto de referencia, los resultados muestran que un bucle for es más eficiente queaccumarray .


Puedes usar findy accumarray:

[c, r] = find(A');
C = accumarray(r, c, [], @(v) {v'});

La matriz se transpone ( A') porque se findagrupa por columna.

Ejemplo:

A = [1 0 0 1 0
     0 1 0 0 0
     0 0 1 1 0
     1 0 1 0 1];

%  Find nonzero rows and colums
[c, r] = find(A');

%  Group row indices for each columns
C = accumarray(r, c, [], @(v) {v'});

% Display cell array contents
celldisp(C)

Salida:

C{1} = 
     1     4

C{2} = 
     2

C{3} =
     3     4

C{4} = 
     1     3     5

Punto de referencia:

m = 10000;
n = 10000;

A = randi([0 1], m,n);

disp('accumarray:')
tic
[c, r] = find(A');
C = accumarray(r, c, [], @(v) {v'});
toc
disp(' ')

disp('For loop:')
tic
C = cell([size(A,1) 1]);
for i = 1:size(A,1)
    C{i} = find(A(i,:));
end
toc

Resultado:

accumarray:
Elapsed time is 2.407773 seconds.

For loop:
Elapsed time is 1.671387 seconds.

Un bucle for es más eficiente que accumarray...

Eliahu Aaron
fuente
Este es más o menos el método ya propuesto por obchardon , ¿no?
Wolfie
Sí, fui un poco lento, vi su respuesta después de publicar la mía.
Eliahu Aaron
2

Usando accumarray :

M = [0 1 0
     1 0 1
     1 1 1];

[val,ind] = find(M.');

A = accumarray(ind,val,[],@(x) {x});
obchardon
fuente
1
El tiempo de ejecución de Octave y MATLAB en línea es de aproximadamente 2 veces de un simple bucle for como: MM{I} = find(M(I, :)).
HansHirse
2
@Hans tal vez quieras ver mi respuesta
Wolfie
Sí, dado que el tamaño de cada celda no es el mismo, este problema no puede ser completamente vectorizado (o hay un truco que no he visto). Es solo una solución que oculta el bucle for.
obchardon
No hay necesidad de ind2sub:[ii, jj] = find(M); accumarray(ii, jj, [], @(x){x})
Luis Mendo
@LuisMendo gracias, he editado mi respuesta.
obchardon
2

Puedes usar strfind :

A = strfind(cellstr(char(M)), char(1));
rahnema1
fuente
Yo (perezosamente) ni siquiera he mirado en los documentos, pero ¿sería más rápido usar stringtipos reales , en lugar de caracteres? Hay muchas optimizaciones para las cadenas, de ahí por qué existen ...
Wolfie
@Wolfie Creo que las matrices numéricas son más similares a las matrices de caracteres que las cadenas, por lo que la conversión de una matriz numérica a una matriz de caracteres debería ser más sencilla que la conversión a una cadena.
rahnema1