Wektoryzacja Matlaba - niezerowe indeksy wierszy macierzy do komórki

10

Pracuję z Matlabem.

Mam binarną macierz kwadratową. Dla każdego wiersza znajduje się jeden lub więcej wpisów 1. Chcę przejrzeć każdy wiersz tej macierzy i zwrócić indeks tych 1s i zapisać je we wpisie komórki.

Zastanawiałem się, czy można to zrobić bez zapętlania wszystkich wierszy tej macierzy, ponieważ w Matlabie pętla jest naprawdę wolna.

Na przykład moja matryca

M = 0 1 0
    1 0 1
    1 1 1 

W końcu chcę coś takiego

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

Podobnie Ajak komórka.

Czy istnieje sposób na osiągnięcie tego celu bez użycia pętli for, w celu szybszego obliczenia wyniku?

ftxx
źródło
Czy chcesz, aby wynik był szybki, czy chcesz uniknąć forpętli? W przypadku tego problemu, w przypadku nowoczesnych wersji MATLAB, podejrzewam, że forpętla będzie najszybszym rozwiązaniem. Jeśli masz problem z wydajnością, podejrzewam, że szukasz niewłaściwego miejsca na rozwiązanie w oparciu o nieaktualne porady.
Czy
@ Chcę, aby wyniki były szybkie. Moja matryca jest bardzo duża. Czas pracy wynosi około 30 sekund na moim komputerze za pomocą pętli for. Chcę wiedzieć, czy są jakieś sprytne operacje wektoryzacji lub mapReduce itp., Które mogą zwiększyć prędkość.
ftxx
1
Podejrzewam, że nie możesz. Wektoryzacja działa na dokładnie opisanych wektorach i macierzach, ale twój wynik pozwala na wektory o różnych długościach. Tak więc, zakładam, że zawsze będziesz mieć jakąś wyraźną pętlę lub coś w rodzaju ukrytej pętli cellfun.
HansHirse
@ftxx jak duży? A ile 1jest w typowym rzędzie? Nie spodziewałbym się, że findpętla zajmie coś blisko 30s dla czegoś wystarczająco małego, aby zmieściło się w pamięci fizycznej.
Czy
@ftxx Zobacz moją zaktualizowaną odpowiedź, którą edytowałem, ponieważ została zaakceptowana z niewielką poprawą wydajności
Wolfie

Odpowiedzi:

11

Na dole tej odpowiedzi znajduje się kod porównawczy, ponieważ wyjaśniłeś, że interesuje Cię wydajność, a nie arbitralne unikanie forpętli.

W rzeczywistości myślę, że forpętle są prawdopodobnie najbardziej wydajną opcją tutaj. Od czasu wprowadzenia „nowego” (2015b) silnika JIT pętle ( źródłowe ) fornie są z natury wolne - w rzeczywistości są zoptymalizowane wewnętrznie.

Można zobaczyć od benchmarku, że mat2cellopcja oferowana przez ThomasIsCoding tutaj jest bardzo powolny ...

Porównanie 1

Jeśli pozbywamy się tej linii, aby skala była wyraźniejsza, wówczas moja splitapplymetoda jest dość powolna, opcja akumulacji tablicy obchardon jest nieco lepsza, ale najszybsze (i porównywalne) opcje albo używają arrayfun(jak sugeruje Thomas), albo forpętli. Pamiętaj, że w większości przypadków arrayfunjest to forukryta pętla, więc nie jest to zaskakujący remis!

Porównanie 2

Polecam użycie forpętli dla zwiększenia czytelności kodu i najlepszej wydajności.

Edytuj :

Jeśli założymy, że zapętlenie jest najszybszym podejściem, możemy dokonać optymalizacji wokół findpolecenia.

konkretnie

  • Zrób Mlogiczne. Jak pokazuje poniższy wykres, może to być szybsze dla stosunkowo małych M, ale wolniejsze z kompromisem konwersji typu dla dużych M.

  • Użyj logiki, Maby zindeksować tablicę 1:size(M,2)zamiast używać find. Pozwala to uniknąć najwolniejszej części pętli ( findpolecenia) i przewyższa narzut związany z konwersją typów, co czyni ją najszybszą opcją.

Oto moja rekomendacja dla najlepszej wydajności:

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

Dodałem to do poniższego testu porównawczego, oto porównanie podejść w stylu pętli:

Porównanie 3

Kod porównawczy:

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
źródło
1
Już widziałem i głosowałem. :-) Wciąż czekam na Luisa; z pewnością ma do tego czarną magię MATLAB.
HansHirse
@Hans Haha tak, chociaż jego zwykły zestaw sztuczek (niejawna ekspansja, sprytne indeksowanie, ...) zwykle utrzymuje rzeczy jako matryce, wąskie gardło tutaj podsumowuje się w komórkach
Wolfie
1
Zauważ, że czasy te są silnie zależne od rzadkości M. Jeśli na przykład zapełnione jest tylko 5% elementów, M = randi([0,20],N) == 20;wówczas forpętla jest zdecydowanie najwolniejsza, a twoja arrayfunmetoda wygrywa.
Czy
@HansHirse :-) Moje podejście byłoby accumarraybez ind2sub, ale jest wolniejsze niż forpętla
Luis Mendo
2

Możesz spróbować arrayfunjak poniżej, które przesuwają się przez rzędyM

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

}

lub (wolniejsze podejście do 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
źródło
1
Chociaż arrayfunjest to po prostu przebranie w pętli, więc może się nie udać na obu frontach 1) unikania pętli i 2) szybkiego, zgodnie z oczekiwaniami OP
Wolfie
2

Edycja : Dodałem test porównawczy, wyniki pokazują, że pętla for jest bardziej wydajna niżaccumarray .


Możesz użyć findi accumarray:

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

Macierz jest transponowana ( A'), ponieważ findgrupuje według kolumn.

Przykład:

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)

Wynik:

C{1} = 
     1     4

C{2} = 
     2

C{3} =
     3     4

C{4} = 
     1     3     5

Reper:

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

Wynik:

accumarray:
Elapsed time is 2.407773 seconds.

For loop:
Elapsed time is 1.671387 seconds.

Pętla for jest bardziej wydajna niż accumarray...

Eliahu Aaron
źródło
Tak, byłem trochę powolny, zobaczyłem jego odpowiedź po opublikowaniu mojej.
Eliahu Aaron
2

Za pomocą accumarray :

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

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

A = accumarray(ind,val,[],@(x) {x});
obchardon
źródło
1
Czas realizacji w oktawie i MATLAB Online jest około 2x prostej pętli for jak: MM{I} = find(M(I, :)).
HansHirse
2
@Hans, możesz chcieć zobaczyć moją odpowiedź
Wolfie
tak, ponieważ rozmiar każdej komórki nie jest taki sam, tego problemu nie można w pełni wektoryzować (lub istnieje sztuczka, której nie widziałem). To tylko rozwiązanie, które ukrywa pętlę for.
obchardon
Nie ma potrzeby ind2sub:[ii, jj] = find(M); accumarray(ii, jj, [], @(x){x})
Luis Mendo
@LuisMendo dzięki, zredagowałem swoją odpowiedź.
obchardon
2

Możesz użyć strfind :

A = strfind(cellstr(char(M)), char(1));
rahnema1
źródło
Ja (leniwie) nawet nie przeglądałem dokumentów, ale czy byłoby to szybsze przy użyciu rzeczywistych stringtypów, a nie znaków? Istnieje wiele optymalizacji dla ciągów, dlatego też istnieją ...
Wolfie,
@Wolfie Myślę, że tablice numeryczne są bardziej podobne do tablic znakowych niż ciągów, więc konwersja tablicy numerycznej na tablicę znaków powinna być prostsza niż konwersja na ciąg znaków.
rahnema1