Pytanie Jak wektoryzować ten kod?


Napisałem funkcję rekursywną, jednak zajmuje ona dużo czasu. Dlatego wektoryzowałem to, ale nie przynosi ono takiego samego rezultatu, jak funkcja rekursywna. To jest mój niezalogowany kod:

function visited = procedure_explore( u, adj_mat, visited )
visited(u) = 1;
neighbours = find(adj_mat(u,:));
for ii = 1:length(neighbours)
    if (visited(neighbours(ii)) == 0)
        visited = procedure_explore( neighbours(ii), adj_mat, visited );
    end
end
end

To jest mój wektoryzowany kod:

function visited = procedure_explore_vec( u, adj_mat, visited )
visited(u) = 1;
neighbours = find(adj_mat(u,:));
len_neighbours=length(neighbours);
visited_neighbours_zero=visited(neighbours(1:len_neighbours)) == 0;
if(~isempty(visited_neighbours_zero))
    visited = procedure_explore_vec( neighbours(visited_neighbours_zero), adj_mat, visited );
end
end

To jest kod testowy

function main
    adj_mat=[0 0 0 0;
             1 0 1 1;
             1 0 0 0;
             1 0 0 1];
    u=2;
    visited=zeros(size(adj_mat,1));
    tic
    visited = procedure_explore( u, adj_mat, visited )
    toc
    visited=zeros(size(adj_mat,1));
    tic
    visited = procedure_explore_vec( u, adj_mat, visited )
    toc
end

Oto algorytm, który próbuję zaimplementować: enter image description here

Jeśli wektoryzacja jest niemożliwa, rozwiązanie mex byłoby również dobre.

Zaktualizuj benchmark: Ten benchmark oparty jest na MATLAB 2017a. Pokazuje, że oryginalny kod jest szybszy niż inne metody

Speed up between original and logical methods is 0.39672
Speed up between original and nearest methods is 0.0042583

Pełny kod

function main_recersive
    adj_mat=[0 0 0 0;
             1 0 1 1;
             1 0 0 0;
             1 0 0 1];
    u=2;
    visited=zeros(size(adj_mat,1));
    f_original=@()(procedure_explore( u, adj_mat, visited ));
    t_original=timeit(f_original);

    f_logical=@()(procedure_explore_logical( u, adj_mat ));
    t_logical=timeit(f_logical);

    f_nearest=@()(procedure_explore_nearest( u, adj_mat,visited ));
    t_nearest=timeit(f_nearest);

    disp(['Speed up between original and logical methods is ',num2str(t_original/t_logical)])
    disp(['Speed up between original and nearest methods is ',num2str(t_original/t_nearest)])    

end

function visited = procedure_explore( u, adj_mat, visited )
    visited(u) = 1;
    neighbours = find(adj_mat(u,:));
    for ii = 1:length(neighbours)
        if (visited(neighbours(ii)) == 0)
            visited = procedure_explore( neighbours(ii), adj_mat, visited );
        end
    end
end

function visited = procedure_explore_nearest( u, adj_mat, visited )
    % add u since your function also includes it.
    nodeIDs = [nearest(digraph(adj_mat),u,inf) ; u];
    % transform to output format of your function
    visited = zeros(size(adj_mat,1));
    visited(nodeIDs) = 1;

end 

function visited = procedure_explore_logical( u, adj_mat )
   visited = false(1, size(adj_mat, 1));
   visited(u) = true;
   new_visited = visited;
   while any(new_visited)
      visited = any([visited; new_visited], 1);
      new_visited = any(adj_mat(new_visited, :), 1);
      new_visited = and(new_visited, ~visited);
   end
end

12
2017-10-04 12:54


pochodzenie


Cóż, twoja funkcja nie jest w stanie działać na wielu wejściach w drugiej wersji, więc przejście w wektorze nie powoduje magicznej wektoryzacji. - Mad Physicist
Ale doskonałe pytanie pod względem klarownego opisu i kompletnego minimalnego przykładu. Spróbuję wymyślić coś na mojej drodze do komputera. - Mad Physicist
Dzięki, Szalony Fizyk. W przypadku gdy wektoryzacja jest niemożliwa, akceptuję również kod wersji mex. - Jame
Nie ma problemu. Zostawię ci mexing, jeśli nic nie wymyślę. W końcu wektoryzacja jest w zasadzie taka sama jak w przypadku mexowania. Po prostu delegujesz pętlę do szybszej implementacji tego samego. - Mad Physicist
Jeśli twoje macierze mają być 4x4, twoje testy są w porządku. Jeśli zamierzasz używać większych wykresów, powinieneś użyć większych macierzy testowych. - beaker


Odpowiedzi:


Oto mała, zabawna funkcja, która przeprowadza nierekursywne wyszukiwanie po raz pierwszy na wykresie.

function visited = procedure_explore_logical( u, adj_mat )
   visited = false(1, size(adj_mat, 1));
   visited(u) = true;
   new_visited = visited;

   while any(new_visited)
      visited = any([visited; new_visited], 1);
      new_visited = any(adj_mat(new_visited, :), 1);
      new_visited = and(new_visited, ~visited);
   end
end

W Octave to działa około 50 razy szybciej niż twoja wersja rekurencyjna na matrycy przyległości 100x100. Musisz przetestować go na MATLAB, aby zobaczyć, co otrzymujesz.


4
2017-10-06 19:34





Możesz myśleć o swojej macierzy sąsiedztwa jako liście ścieżek o długości dokładnie jeden. Możesz generować ścieżki o innych długościach n, przenosząc je do n-tej potęgi do rangi macierzy. (adj_mat ^ 0 to macierz tożsamości)

Na wykresie z n węzłów najdłuższa ścieżka nie może być dłuższa niż n-1, dlatego można zsumować wszystkie moce dla analizy osiągalności:

adj_mat + adj_mat^2 + adj_mat^3
ans =
   0   0   0   0
   4   0   1   3
   1   0   0   0
   3   0   0   3

Jest to liczba (różnych) sposobów, za pomocą których można przejść z jednego węzła do drugiego. Aby uzyskać prosty zasięg, sprawdź, czy ta wartość jest większa od zera:

visited(v) = ans(v, :) > 0;

W zależności od definicji może zajść konieczność zmiany kolumn i wierszy w wyniku (np. Take ans (:, v)).

Aby zwiększyć wydajność, możesz użyć niższych mocy, aby uzyskać wyższe. Na przykład coś takiego jak A + A ^ 2 + A ^ 3 + A ^ 4 + A ^ 5 byłoby wydajnie obliczone:

A2 = A^2;
A3 = A2*A
A4 = A2^2;
A5 = A4*A;
allWalks= A + A2 + A3 + A4 + A5;

Uwaga: Jeśli chcesz dołączyć oryginalny węzeł jako osiągalny z samego siebie, powinieneś dołączyć do sumy macierz tożsamości.

Minimalizuje to liczbę multiplikacji macierzy, a także MATLAB prawdopodobnie wykona kwadraty macierzy szybciej niż regularne mnożenie.

Z mojego doświadczenia wynika, że ​​mnożenie macierzy jest stosunkowo szybkie w MATLAB-ie, co da wynik równy wynikowi (osiągalności) dla wszystkich węzłów na wykresie naraz. Jeśli interesuje Cię tylko niewielki podzbiór dużego wykresu, prawdopodobnie nie jest to najlepsze rozwiązanie.

Zobacz także tę odpowiedź: https://stackoverflow.com/a/7276595/1974021


2
2017-10-06 20:51



Cześć, co jest v? - Jame
Wziąłem v z twojego pytania "odwiedzony (u) jest ustawiony na true dla wszystkich węzłów dostępnych z v" - DasKrümelmonster


Nie sądzę, żebyś mógł odpowiednio wektoryzować swoją funkcję: Twoja oryginalna funkcja nigdy nie dociera wielokrotnie do tego samego węzła. Wektoryzując go, przeszedłbyś wszystkie bezpośrednio połączone węzły w tym samym czasie do następnej funkcji. Dlatego możliwe byłoby, że w kolejnych przypadkach ten sam węzeł zostanie osiągnięty wiele razy. Na przykład. w twoim przykładzie węzeł 1 zostanie osiągnięty 3 razy. Tak więc, chociaż nie ma już pętli, funkcja może, w zależności od sieci, być wywoływana więcej razy rekurencyjnie, co może wydłużyć czas obliczeniowy.

Biorąc to pod uwagę, nie jest ogólnie możliwe znalezienie wszystkich dostępnych węzłów bez pętli lub wywołań rekursywnych. Na przykład możesz sprawdzić wszystko (poprawne lub nieprawidłowe) ścieżki. Ale to działałoby zupełnie inaczej niż twoja funkcja i, w zależności od liczby węzłów, może spowodować utratę wydajności ze względu na oszałamiającą liczbę ścieżek do sprawdzenia. Twoja obecna funkcja nie jest zbyt zła i będzie się dobrze skalować w dużych sieciach.

Trochę offtopic, ale od Matlab 2016a możesz go używać nearest() aby znaleźć wszystkie osiągalne węzły (bez początkowego węzła). Wywołuje algorytm pierwszy w stosunku do pierwszego algorytmu głębi:

% add u since your function also includes it.
nodeIDs = [nearest(digraph(adj_mat),u,inf) ; u]; 

% transform to output format of your function
visited = zeros(size(adj_mat,1));
visited(nodeIDs) = 1;

Jeśli jest to projekt dla studentów, możesz argumentować, że podczas gdy twoja funkcja działa, użyłeś wbudowanej funkcji ze względu na wydajność.


1
2017-10-05 10:04



Przepraszamy za odpowiedź z opóźnieniem. Muszę zainstalować Malab 2017b do testowania twojego kodu. Udało się, ale czasochłonny jest większy niż mój oryginalny kod. To złe wieści. Jednak dzięki za rozwiązanie. Nie znałem wbudowanej funkcji. Będę czekać na inne rozwiązanie (mex), które może być szybsze niż rozwiązanie i moje rozwiązanie - Jame
Wow, możesz być dumny, że twoja funkcja jest szybsza niż funkcja wbudowana - to raczej rzadkość! Zawsze zakładam, że jakakolwiek funkcja chłopaków z Mathworks jest o wiele rzędów wielkości szybsza od czegoś napisanego od zera. - Leander Moesinger
Tak. Możesz to sprawdzić. Użyłem Matlab 2017b :) - Jame
@Jame Jeśli osiągnąłeś ten rodzaj wykonania, zazwyczaj można bezpiecznie powiedzieć, że zmiana kodu nie spowoduje znaczących ulepszeń. Tutaj może nie być to możliwe, ale w tych przypadkach wdrożenie innego rozwiązania (/ przyjrzenie się problemowi root) jest zwykle jedynym sposobem na prawdziwy krok naprzód. - Dennis Jaheruddin


Problem z funkcją rekurencyjną jest związany z visited(u) = 1;. Dzieje się tak, ponieważ MATLAB wykorzystuje technikę kopiowania przy zapisie, aby przekazywać / przypisać zmienne. Jeśli się nie zmienisz visited w treści funkcji nie jest wykonywana żadna kopia, ale gdy jest modyfikowana, tworzona jest jej kopia, a modyfikacja jest stosowana na jej kopii. Aby temu zapobiec, możesz użyć a uchwyt obiektu do przekazania przez odniesienie do funkcji.

Zdefiniuj klasę uchwytu (zapisz ją na visited_class.m):

classdef visited_class < handle
    properties
        visited
    end
    methods
        function obj = visited_class(adj_mat)
            obj.visited = zeros(1, size(adj_mat,1));
        end
    end
end

Funkcja rekursywna:

function procedure_explore_handle( u, adj_mat,visited_handle )
    visited_handle.visited(u) = 1;
    neighbours = find(adj_mat(u,:));
    for n = neighbours
        if (visited_handle.visited(n) == 0)
            procedure_explore_handle( n, adj_mat , visited_handle );
        end
    end
end

Inicjalizuj zmienne:

adj_mat=[0 0 0 0;
         1 0 1 1;
         1 0 0 0;
         1 0 0 1];
visited_handle = visited_class(adj_mat);
u = 2;

Nazwij to jako:

procedure_explore_handle( u, adj_mat,visited_handle );

Wynik zostanie zapisany w visited_handle:

disp(visited_handle.visited)

1
2017-10-12 07:15





Jeśli chcesz przejść od jednego punktu wykresu do drugiego, najskuteczniejszym sposobem na znalezienie go pod względem zasobów jest algorytm Dijkstry. Algorytm Floyd-Warshall oblicza wszystkie odległości między wszystkimi punktami i może być równoległy (rozpoczynając od wielu punktów).

Dlaczego trzeba wektoryzować (lub używać programowania mex)? Jeśli chcesz po prostu w jak największym stopniu korzystać z szybkich procesów mnożenia macierzy Matlaba, użycie produktów A powinno Cię tam szybko:

adj_mat2=adj_mat^2;               % allowed to use 2 steps
while (adj_mat2 ~= adj_mat)       % check if new points were reached
      adj_mat=adj_mat2;           % current set of reachable points
      adj_mat2=(adj_mat^2)>0;     % allowed all steps again: power method
end

0
2017-10-09 08:29



Gdzie skorzystałeś u? - Jame
O przepraszam. Jeśli macierz sąsiedztwa jest symetryczna u=(sum(adj_mat2))>0 po zakończeniu pętli. - alle_meije


Ta odpowiedź daje tylko jednoznaczną, wektoryzowaną implementację sugestii od Odpowiedź DasKrümelmonster, które moim zdaniem jest szybsze niż kod w pytaniu (przynajmniej jeśli wymiary macierzy nie są zbyt duże). Korzysta z polyvalm funkcja do oceny sumy mocy macierzy sąsiedztwa.

function visited = procedure_explore_vec(u, adj_mat)
    connectivity_matrix = polyvalm(ones(size(adj_mat,1),1),adj_mat)>0;

    visited = connectivity_matrix(u,:);
end

0
2017-10-11 17:25