Wyszukiwanie maksimum lokalnego cz II. - O szybkim turyście i chciwym bacy

Treść tego zadania jest bardzo głupia :) Chciałem ją zmodyfikować, ale stwierdziłem, że tak będzie śmieszniej :) Miłej lektury.

Mamy dany graf w którym wyróżniono dwa wierzchołki - A oraz B. Z każdą krawędzią i są związane dwie liczby rzeczywiste: a_i oraz b_i. Koszt przejechania tej krawędzi w czasie t wynosi a_it + b_i.
Turysta chce przejść z punktu A do B. Turysta jest bardzo sprytny i w dowolnym czasie t potrafi znaleźć najtańszą drogę z punktu A do B. Jest też bardzo bardzo szybki (struś pędziwiatr to przy nim żółw) i całą trasę potrafi przebiec w czasie zerowym!

Pytanie: W jakim czasie T z przedziału [0, 1] baca powinien wypuścić turystę ze swojej chatki, aby zmaksymalizować swój zysk?

Co ma z tym wspólnego maksimum lokalne? Dla każdej ścieżki w grafie jej koszt jest funkcją liniową zależną od T.
Zatem baca jeśli wypuści turystę w czasie T zarobi min(T) dudków, gdzie min(T) jest dane wzorem:
Spójrzmy na wykres tej funkcji.
Funkcja ma tylko jedno maksimum lokalne. Możemy zatem użyć wyszukiwania ternarnego (wystarczy tylko uogólnić algorytm na przedział rzeczywisty - co pozostawiam jako proste ćwiczenie ;) )

Wyszukiwanie maksimum lokalnego cz I. - wyszukiwanie ternarne

Mamy daną tablicę A[1..n] wypełnioną parami różnymi liczbami rzeczywistymi. Element, który jest większy od swoich sąsiadów nazywamy maksimum lokalnym. W tej notce przedstawię algorytm, który dla danej tablicy znajdzie (jakieś) maksimum lokalne.
W pierwszej kolejności sprawdźmy elementy A[1] oraz A[n]. Mają one tylko jednego sąsiada. Jeśli zatem A[1] jest większe niż A[2] albo A[n] jest większe od A[n-1] to znaleźliśmy maksimum lokalne i kończymy algorytm.
Co jednak jeśli tak nie jest? Wyobraźmy sobie tablice jako wierzchołki. Krawędź między dwoma wierzchołkami jest skierowana w stronę większego elementu. Mamy zatem taką sytuację:
Algorytm działa rekurencyjnie. Oznaczmy zatem element najbardziej na lewo jako A[l], a najbardziej na prawo jako A[r]. Niezmiennikiem naszego wywołania rekurencyjnego będzie to że element A[l] jest mniejszy od A[l+1] a element A[r] jest mniejszy od elementu A[r-1] (rysunek powyżej). Można łatwo udowodnić że w takim przedziale znajduje się maksimum lokalne.
Weźmy dwa środkowe elementy tablicy. Oznaczmy je jako A[p] oraz A[q]. Zapytajmy się który z nich jest większy.
Jeśli A[p] < A[q] to rozwiązanie znajduje się w przedziale A[p..r].
Jeśli A[p] > A[q] to rozwiązanie znajduje się w przedziale A[l..q].
Nie zdarzy się, ża A[p] = A[q] bo założyliśmy że elementy są parami różne. Algorytm zatrzymuje się gdy w dostajemy tablicę 3 elementową (można łatwo udowodnić, że na taką tablicę trafimy: gdy natrafiamy na tablicę rozmiaru 4 to następnie uruchamiamy się dla tablicy 3 elementowej; gdy natrafimy na tablicę rozmiaru 5 to następnie uruchamiamy się dla tablicy albo 3 elementowej albo 4 elementowej). Gdy natrafimy na tablicę 3 elementową to zgodnie z niezmiennikiem algorytmu jego środkowy element będzie maksimum lokalnym.
Algorytm jako że w każdej iteracji wykonuje stałą liczbę operacji, a zawsze przedział dzielony jest na "prawie" połowę to działa w czasie O(lg n).
Jeśli elementy w tablicy najpierw rosną, a potem maleją (lub odwrotnie) to tablica taka ma tylko jedno maksimum (minimum) lokalne i jest nim największy (najmniejszy) element w tablicy. Zatem w tego rodzajach tablic jesteśmy w stanie znaleźć element największy (najmniejszy) w czasie logarytmicznym!

Algorytm Forda-Johnsona cz III.

Zgodnie z obietnicą w tej części postaram się wyjaśnić dlaczego akurat taki ciąg został wybrany w części drugiej.

Wszystkie rozważenia przeprowadzę na następującym przykładzie:Rozważymy trzy proste ciągi i policzymy ile w najgorszym przypadku będziemy potrzebowali porównań, aby posortować te liczby. Przypomnę jeszcze, że aby do ciągu n elementowego wstawić liczbę potrzebujemy [lg n] + 1 porównań (gdzie [x] - oznacza podłogę z liczby x).

Spróbujmy wstawiać elementy od lewej do prawej:

Wstawiamy 8 --- Wykonujemy 4 porównania ([lg 8] + 1)
Wstawiamy 7 --- Wykonujemy 4 porównania ([lg 8] + 1)
Wstawiamy 6 --- Wykonujemy 4 porównania ([lg 8] + 1)
Wstawiamy 5 --- Wykonujemy 4 porównania ([lg 8] + 1)
Wstawiamy 4 --- Wykonujemy 4 porównania ([lg 8] + 1)
Wstawiamy 3 --- Wykonujemy 4 porównania ([lg 8] + 1)
Wstawiamy 2 --- Wykonujemy 4 porównania ([lg 8] + 1)
W sumie - 28 porównań.

Zauważmy, że w tej metodzie wstawiając każdy element, wykonaliśmy tyle samo porównań (skorzystamy z tego faktu). Teraz spróbujmy wstawiać elementy od prawej do lewej:

Wstawiamy 2 --- Wykonujemy 2 porównania ([lg 2] + 1)
Wstawiamy 3 --- Wykonujemy 3 porównania ([lg 4] + 1)
Wstawiamy 4 --- Wykonujemy 3 porównania ([lg 6] + 1)
Wstawiamy 5 --- Wykonujemy 4 porównania ([lg 8] + 1)
Wstawiamy 6 --- Wykonujemy 4 porównania ([lg 10] + 1)
Wstawiamy 7 --- Wykonujemy 4 porównania ([lg 12] + 1)
Wstawiamy 8 --- Wykonujemy 4 porównania ([lg 14] + 1)
W sumie - 24 porównania.

Czyli w drugim przypadku wykonaliśmy o 4 porównania mniej. I nasza metoda z poprzedniej notki:

Wstawiamy 3 --- Wykonujemy 2 porównania ([lg 3] + 1)
Wstawiamy 2 --- Wykonujemy 2 porównania ([lg 3] + 1)
Wstawiamy 5 --- Wykonujemy 3 porównania ([lg 7] + 1)
Wstawiamy 4 --- Wykonujemy 3 porównania ([lg 7] + 1)
Wstawiamy 8 --- Wykonujemy 4 porównania ([lg 12] + 1)
Wstawiamy 7 --- Wykonujemy 4 porównania ([lg 12] + 1)
Wstawiamy 6 --- Wykonujemy 4 porównania ([lg 12] + 1)
W sumie - 22 porównania.

Zauważmy dwie rzeczy. W momentach w którym nasz ciąg maleje - (np. we wstawieniu 8, 7, 6) wykonujemy tyle samo porównań (zauważyliśmy już to w ciągu pierwszym). Nie przypadkiem jest również to, że wewnątrz logarytmów są liczby bliskie potęgom liczby 2 (gdybyśmy wstawiali liczby do 11 to w miejscu [lg 12] pojawił by się [lg 15]).

Ciąg trzeci został więc wybrany następująco (zachłannie):
- wybierzmy wierzchołek najbardziej na lewo, który jesteśmy w stanie wstawić w dwóch ruchach.
- wstawmy go do ciągu "głównego" wraz z wszystkimi elementami na prawo od niego
- wybierzmy wierzchołek najbardziej na lewo, który jesteśmy w stanie wstawić w trzech ruchach.
- wstawmy go do ciągu "głównego" wraz z wszystkimi elementami na prawo od niego
- wybierzmy wierzchołek najbardziej na lewo, który jesteśmy w stanie wstawić w czterech ruchach.
- wstawmy go do ciągu "głównego" wraz z wszystkimi elementami na prawo od niego
...

Mam nadzieję, że nieco się wyjaśniło odnośnie tego ciągu. Jeśli nie to proszę o kontakt :)

Algorytm Forda-Johnsona cz II.

Algorytm Forda Johnsona służy do posortowania n-elementowej tablicy. Algorytm jest bardzo ciekawy z teoretycznego punktu widzenia, bo wykonuje bardzo małą liczbę porównań, bliską optymalnej. Algorytm działa optymalnie dla n < 16 oraz dla n = 20, 21, 22, a najmniejszy n dla którego znamy algorytm wymagający mniejszej liczby porównań to n = 47.

Na początku podzielmy tablicę n-elementową na pary i porównujemy je w tych parach. Jeśli liczba elementów jest nieparzysta to jeden element zostaje nieporównany.Elementy, które okazały się mniejsze w parach, sortujemy rekurencyjnie. Dla ułatwienia naszych dalszych rozważań ponumerujmy elementy. Największy numer dostaje element, który nie został jeszcze porównany (jeśli tablica ma nieparzyście wiele elementów). Resztę elementów numerujemy tak jak na rysunku.Teraz każdy element znajdujący się w ciągu "na dole" będziemy chcieli umieścić w ciągu "na górze". Szukając miejsca tego elementu, będziemy korzystać z wyszukiwania binarnego.
Kolejność w której będziemy wrzucać elementy z dołu do elementów na górze ma kolosalne znaczenie. Optymalną kolejnością wrzucania elementów jest następująca: 1, 3, 2, 5, 4, 11, 10, 9, 8, 7, 6, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 43, 42 ... . Element 1 wpisałem dla porządku, element ten jest już na swoim miejscu i nie trzeba nic z nim robić. Ciąg niemniej może wydawać się trochę skomplikowany. Już tłumaczę jak powstał. Można zauważyć, że po wszystkich liczbach zaznaczonych pogrubieniem, liczby maleją, aż skończą się nam liczby. Następnie pojawia się znowu pogrubiona liczba i liczby znowu zmniejszają się o jeden. Aby umieć wygenerować taki ciąg wystarczy zatem wiedzieć na jakiej podstawie wybieramy pogrubione liczby. Zauważmy, że sumy kolejnych dwóch liczb pogrubionych są równe kolejnym potęgom dwójki: 1 + 3 = 4; 3 + 5 = 8; 5 + 11 = 16; 11 + 21 = 32 itd.

Algorytm jest nieco skomplikowany. Jeśli są z nim jakieś problemy, albo trzeba coś jeszcze wytłumaczyć - to proszę o kontakt.

Algorytm Forda-Johnsona cz I.

Spróbujmy posortować ciąg 5 liczb za pomocą 7 porównań. Oznaczmy te elementy jako a, b, c, d oraz e. Porównajmy element a z elementem b i załóżmy, że a jest mniejszy (ew. zamieńmy a z b miejscami). Porównajmy element c z d i przyjmijmy, że element c jest mniejszy. I wreszcie porównajmy element a z c i przyjmijmy, że element a jest mniejszy. Wykonaliśmy do tej pory 3 porównania, a sytuację w której się znaleźliśmy przedstawia poniższy rysunek.Spróbujmy teraz wcisnąć element e gdzieś do ciągu [a, c, d]. Możemy to zrobić za pomocą dwóch porównań (użyjemy wyszukiwania binarnego). Najpierw element e sprawdzamy z elementem c. Jeśli jest mniejszy to sprawdzamy z elementem a; w przeciwnym przypadku z elementem d. Trudno przewidzieć w którym miejscu stanie element e. Rozważmy dwa przypadki. Albo e jest najmniejszym elementem z tej trójki (a, c, d), albo nie jest. Ten drugi przypadek przedstawia poniższy rysunek.Jedyną rzeczą jaką trzeba teraz zrobić to wstawić element b pomiędzy elementy {c, d, e}. Identycznie jak poprzednio możemy to zrobić w dwóch ruchach.
W drugim przypadku, gdy element e jest najmniejszy, wystarczy element b porównać tylko z elementami c oraz d. To też możemy zrobić używając dwóch ruchów.

Reasumując. Na początku porównując elementy a, b, c oraz d - wykonaliśmy 3 porównania. Wpychając element e pomiędzy elementy [a, c, e] wykonaliśmy dodatkowe 2 porównania. Na końcu ustawiliśmy element b na swoje miejsce używając nie więcej niż 2 porównań. 3 + 2 + 2 = 7

W następnej notce uogólnimy ten algorytm na tablicę n-elementową.

Wyszukiwanie interpolacyjne

Mamy daną posortowaną N-elementową tablicę zawierającą losowe liczby. Chcemy sprawdzić czy znajduje się wśród nich liczba k. Niejeden programista w tym miejscu krzyknie: "Wyszukiwanie binarne!". Pytanie: Czy można szybciej? Można!

Przypomnijmy może szybko implementację wyszukiwania binarnego:

int bsearch(int* A, int N, int k)
{
  int l = 0, r = N-1;

  while (l <= r)
  {
    int m = (l + r) / 2;

    if (A[m] == k) return m;
    else if(A[m] < k) l = m + 1;
    else r = m - 1;
  }

  return -1;
}

W poszukiwaniu interpolacyjnym korzystamy z następującego spostrzeżenia. Gdy przeszukujemy książkę telefoniczną szukając pana Wołodyjowskiego, to zazwyczaj nie otwieramy książki na stronie środkowej (jak to robi wyszukiwanie binarne) lecz na jej końcu. Z kolei gdy szukamy w książce pani Adamskiej to zaglądamy na początek książki. Aby zaimplementować ten algorytm, wystarczy w algorytmie wyszukiwania binarnego zamienić linijkę:

int m = (l + r) / 2;

na następującą:

int m = l + (k - A[l]) * (r - l) / (A[r] - A[l]);

Zaletą takiego rozwiązania jest złożoność. Algorytm dla tablicy wypełnionej losowymi wartościami używa mniej niż lg lg N + 1 porównań (dla przypomnienia - wyszukiwanie binarne wykonuje lg N + 1 porównań). Funkcja ta rośnie bardzo wolno. Dla przykładu lg lg 10³⁰⁸ < 10. A w przypadku, gdy tablica jest jeszcze regularniej ułożona niż w sposób losowy to liczba wymaganych porównań jest jeszcze mniejsza!

Wadą takiego rozwiązania jest koszt obliczenia wartości m (wykonujemy tutaj dużo operacji arytmetycznych) co dla nie wielkich wartości N może spowodować, że algorytm wyszukiwania binarnego będzie działał szybciej.

Ponadto algorytm nie nadaje się do wyszukiwania w tablicach, o których nie wiemy czy elementy w niej zawarte są równomiernie rozłożone. Algorytm dla przypadków pesymistycznych może wymagać nawet czasu liniowego.

Wzory Cramera

Powiedzmy, że mamy dany układ 2 (3) równań liniowych z 2 (3) niewiadomymi o współczynnikach całkowitych. Wiemy, że układ ma dokładnie jedno rozwiązanie. Chcemy je obliczyć, przy czym chcemy liczyć na liczbach całkowitych tak długo jak tylko się da. Świetnie do tego nadają się wzory Cramera.
Niech dany będzie następujący układ równań:
Macierzą A nazywamy macierzą współczynników tego układu, a macierz B macierzą wyrazów wolnych.
Macierz G_i powstaje w ten sposób, że w macierzy A kolumnę i-tą zastępujemy macierzą B. Dla przykładu macierz G₂ będzie wyglądał następująco:
Rozwiązanie naszego równania wygląda teraz następująco:
Gdzie det X oznacza wyznacznik macierzy X (http://pl.wikipedia.org/wiki/Wyznacznik_macierzy). Jeżeli wyznacznik macierzy A jest równy zero to układ równań nie ma rozwiązania (lub ma ich nieskończenie wiele).
Algorytm wyznacznika macierzy nie zawsze jest łatwa do implementacji. Jednak dla macierzy 2x2 oraz 3x3 obliczenie wyznacznika nie jest trudnym zadaniem.
Co ważne: podczas obliczania wyznacznika korzystamy z operacji dodawania, odejmowania i mnożenia. Nie wychodzi zatem poza liczby całkowite. Operację dzielenia dokonujemy dopiero na samym końcu.
Na koniec jeszcze powiem, że wzory Cramera działają dla każdego układu n równań z n niewiadomymi. Założenia o tym, że mamy tylko 2 (3) równania oraz że współczynniki są całkowite, poczyniliśmy aby pokazać w jakich sytuacjach najlepiej korzystać z metody tu przedstawionej.

Wierszami czy kolumnami? cz III. (Mnożenie macierzy)

W niniejszej notce omówię algorytm mnożenia macierzy. Jeśli czytelnik z mnożeniem macierzy spotyka się po raz pierwszy to odsyłam do wikipedii (http://pl.wikipedia.org/wiki/Mno%C5%BCenie_macierzy). Operacja mnożenia macierzy jest tam bardzo ładnie wytłumaczona.
Przed napisaniem tej notki wyszukałem w googlach hasło "mnożenie macierzy C++". Pobieżne przejrzenie wyników potwierdziło moje przypuszczenia - zdecydowana większość programistów pisze pętle w złej kolejności. Przedstawmy najpierw typową implementację mnożenia macierzy:

const int N = 1000;
static int T[N][N];
static int A[N][N];
static int B[N][N];

for (int i = 0; i < N; i++)
  for (int j = 0; j < N; j++)
  {
    T[i][j] = 0;
    for (int k = 0; k < N; k++)
      T[i][j] += A[i][k] * B[k][j];
  }

Gdzie leży błąd można na pierwszy rzut oka nie zauważyć. Wytłumaczę to zatem, aby wszyscy mogli to zauważyć, a w razie innego algorytmu - powtórzyć rozumowanie.
Z tablicą T nie ma problemów. Wewnątrz pętli odwołujemy się do elementu T[i][j]. Przy czym przechodzimy ją po wierszach. To znaczy, że pętla iterująca po j-tach znajduje się wewnątrz pętli iterującej po zmiennej i. Podobna sytuacja jest z tablicą A. Odwołujemy się do elementu A[i][k], czyli przechodzimy ją po wierszach (pętla iterująca po zmiennej k znajduje się wewnątrz pętli iterującej po zmiennej i). Problem jest z tablicą B.
W przypadku tablicy B odwołujemy się do elementu B[k][j]. Przechodzimy ją zatem po kolumnach - pętla iterująca po j-tach znajduje się na zewnątrz pętli iterującej po zmiennej k. Procesor zatem nie będzie wczytywał optymalnie tablicy B do pamięci cache.
Możemy to poprawić przestawiając miejscami pętle ze zmienną k z pętlą ze zmienną j. Musimy też stworzyć osobne pętle zerujące tablicę.

const int N = 1000;
static int T[N][N];
static int A[N][N];
static int B[N][N];

for (int i = 0; i < N; i++)
  for (int j = 0; j < N; j++)
    T[i][j] = 0;

for (int i = 0; i < N; i++)
  for (int k = 0; k < N; k++)
    for (int j = 0; j < N; j++)
      T[i][j] += A[i][k] * B[k][j];

Tutaj już wszystko powinno działać w porządku. Wszystkie tablice przechodzimy po wierszach.
Przykład ten jest często wykorzystywany przy omawianiu przechodzenia tablicy po wierszach i kolumnach. Oba algorytmy mają nawet swoje nazwy potoczne. Pierwszy z nich nazywamy algorytmem ijk, natomiast drugi (jak nietrudno się domyślić) ikj.
Czas na przetestowanie obu rozwiązań. Różnica jest odczuwalna:
algorytm ijk - 0:13.63
algorytm ikj - 0:10.45
Motto tego trzy częściowego kursu jest następujące: "Iterowanie po wierszach przyśpiesza program". Czasem o kilka sekund, czasem nawet ośmiokrotnie. Warto o tym pamiętać bo przestawianie dwóch pętli nas nic nie kosztuje.