Vectorisation Matlab - indices de ligne de matrice non nulles à la cellule

10

Je travaille avec Matlab.

J'ai une matrice carrée binaire. Pour chaque ligne, il y a une ou plusieurs entrées de 1. Je veux parcourir chaque ligne de cette matrice et retourner l'index de ces 1 et les stocker dans l'entrée d'une cellule.

Je me demandais s'il y avait un moyen de le faire sans boucler sur toutes les lignes de cette matrice, car la boucle est vraiment lente dans Matlab.

Par exemple, ma matrice

M = 0 1 0
    1 0 1
    1 1 1 

Finalement, je veux quelque chose comme

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

Il en Ava de même pour une cellule.

Existe-t-il un moyen d'atteindre cet objectif sans utiliser de boucle for, dans le but de calculer le résultat plus rapidement?

ftxx
la source
Voulez-vous que le résultat soit rapide ou voulez-vous que le résultat évite les forboucles? Pour ce problème, avec les versions modernes de MATLAB, je soupçonne fortement qu'une forboucle est la solution la plus rapide. Si vous avez un problème de performances, je soupçonne que vous cherchez au mauvais endroit la solution basée sur des conseils obsolètes.
Will
@Voudrai-je que les résultats soient rapides. Ma matrice est très grande. Le temps d'exécution est d'environ 30 secondes sur mon ordinateur en utilisant la boucle for. Je veux savoir s'il existe des opérations de vectorisation intelligentes ou mapReduce, etc. qui peuvent augmenter la vitesse.
ftxx
1
Je suppose que tu ne peux pas. La vectorisation fonctionne sur des vecteurs et matrices décrits avec précision, mais votre résultat autorise des vecteurs de longueurs différentes. Ainsi, mon hypothèse est que vous aurez toujours une boucle explicite ou une boucle déguisée comme cellfun.
HansHirse
@ftxx quelle taille? Et combien de 1s dans une rangée typique? Je ne m'attendrais pas à ce qu'une findboucle prenne quelque chose près de 30 secondes pour quelque chose d'assez petit pour tenir sur la mémoire physique.
Le
@ftxx S'il vous plaît voir ma réponse mise à jour, je l'ai modifiée depuis qu'elle a été acceptée avec une amélioration mineure des performances
Wolfie

Réponses:

11

Au bas de cette réponse se trouve un code de référence, car vous avez précisé que vous êtes intéressé par les performances plutôt que d'éviter arbitrairement les forboucles.

En fait, je pense que les forboucles sont probablement l'option la plus performante ici. Depuis que le "nouveau" moteur JIT (2015b) a été introduit ( source ), les forboucles ne sont pas intrinsèquement lentes - en fait, elles sont optimisées en interne.

Vous pouvez voir dans le benchmark que l' mat2celloption offerte par ThomasIsCoding ici est très lente ...

Comparaison 1

Si nous nous débarrassons de cette ligne pour rendre l'échelle plus claire, alors ma splitapplyméthode est assez lente, l' option accumarray d' Obchardon est un peu meilleure, mais les options les plus rapides (et comparables) utilisent soit arrayfun(comme l'a également suggéré Thomas) soit une forboucle. Notez que arrayfunc'est essentiellement une forboucle déguisée pour la plupart des cas d'utilisation, donc ce n'est pas une cravate surprenante!

Comparaison 2

Je vous recommande d'utiliser une forboucle pour une meilleure lisibilité du code et les meilleures performances.

Modifier :

Si nous supposons que le bouclage est l'approche la plus rapide, nous pouvons faire quelques optimisations autour de la findcommande.

Plus précisément

  • Rendez Mlogique. Comme le montre le graphique ci-dessous, cela peut être plus rapide pour les relativement petits M, mais plus lent avec le compromis de conversion de type pour les grands M.

  • Utilisez une logique Mpour indexer un tableau 1:size(M,2)au lieu d'utiliser find. Cela évite la partie la plus lente de la boucle (la findcommande) et l'emporte sur la surcharge de conversion de type, ce qui en fait l'option la plus rapide.

Voici ma recommandation pour de meilleures performances:

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

Je l'ai ajouté au benchmark ci-dessous, voici la comparaison des approches de style boucle:

Comparaison 3

Code de référence:

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
la source
1
Déjà vu et voté. :-) Toujours en attente de Luis; il a certainement de la magie noire MATLAB pour ça.
HansHirse
@Hans Haha oui, bien que son sac de trucs habituel (expansion implicite, indexation intelligente, ...) conserve généralement les choses sous forme de matrices, le goulot d'étranglement ici se résume dans les cellules
Wolfie
1
Notez que ces temps dépendent fortement de la rareté de M. Si, par exemple, seulement 5% des éléments sont remplis, M = randi([0,20],N) == 20;la forboucle est de loin la plus lente et votre arrayfunméthode gagne.
Le
@HansHirse :-) Mon approche aurait été accumarraysans ind2sub, mais elle est plus lente que la forboucle
Luis Mendo
2

Vous pouvez essayer arrayfuncomme ci-dessous, qui parcourent des rangées 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

}

ou (une approche plus lente de 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
la source
1
Bien qu'il s'agisse arrayfunessentiellement d'un déguisement en boucle, cela peut échouer sur les deux fronts: 1) éviter les boucles et 2) être rapide, comme l'espérait l'OP
Wolfie
2

Edit : j'ai ajouté un benchmark, les résultats montrent qu'une boucle for est plus efficace queaccumarray .


Vous pouvez utiliser findet accumarray:

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

La matrice est transposée ( A') car findregroupée par colonne.

Exemple:

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)

Production:

C{1} = 
     1     4

C{2} = 
     2

C{3} =
     3     4

C{4} = 
     1     3     5

Référence:

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

Résultat:

accumarray:
Elapsed time is 2.407773 seconds.

For loop:
Elapsed time is 1.671387 seconds.

Une boucle for est plus efficace que accumarray...

Eliahu Aaron
la source
C'est à peu près la méthode déjà proposée par obchardon , non?
Wolfie
Oui, j'étais un peu lent, j'ai vu sa réponse après avoir posté la mienne.
Eliahu Aaron
2

Utilisation de cumularray :

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

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

A = accumarray(ind,val,[],@(x) {x});
obchardon
la source
1
Temps d'exécution dans Octave et Matlab en ligne est d' environ 2x d'une simple boucle comme: MM{I} = find(M(I, :)).
HansHirse
2
@Hans vous voudrez peut-être voir ma réponse
Wolfie
oui, puisque la taille de chaque cellule n'est pas la même, ce problème ne peut pas être entièrement vectorisé (ou il y a une astuce que je n'ai pas vue). Ce n'est qu'une solution qui cache la boucle for.
obchardon
Pas besoin de ind2sub:[ii, jj] = find(M); accumarray(ii, jj, [], @(x){x})
Luis Mendo
@LuisMendo merci, j'ai édité ma réponse.
obchardon
2

Vous pouvez utiliser strfind :

A = strfind(cellstr(char(M)), char(1));
rahnema1
la source
Je (paresseusement) n'ai même pas regardé dans les documents, mais serait-ce plus rapide d'utiliser des stringtypes réels plutôt que des caractères? Il y a beaucoup d'optimisations pour les chaînes, d'où leur existence ...
Wolfie
@Wolfie Je pense que les tableaux numériques sont plus similaires aux tableaux de caractères qu'aux chaînes, donc la conversion du tableau numérique en tableau de caractères devrait être plus simple que la conversion en chaîne.
rahnema1