O wektoryzacji raz jeszcze

https://xpil.eu/gNVza

Miała być zagadka, ale zamiast tego zrobił mi się temat na Pchełkę, więc dziś pokażę raz jeszcze zalety wektoryzacji.

Mamy do rozwiązania następujące zadanie: bierzemy odcinek o długości 1 i dzielimy go na dwie części w przypadkowym miejscu, z rozkładem równomiernym. Długości obu odcinków mnożymy jedna przez drugą. Następnie powtarzamy tę samą operację dla każdego z tych dwóch odcinków (a więc dzielimy każdy z nich na dwa odcinki w losowym miejscu i mnożymy - osobno dla każdego z nich - długości w każdej parze). Powstałe w ten sposób cztery odcinki poddajemy tej samej procedurze. Powstałe w ten sposób osiem odcinków... i tak dalej.

Po wykonaniu nieskończenie wielu takich podziałów i mnożeń, wszystkie otrzymane iloczyny sumujemy. Pytanie: ile wynosi ta suma?

Z zadaniem można się pewnie zmierzyć algebraicznie (jeżeli są tu jacyś odważni, zapraszam), ja zaatakowałem je symulacją bo na podejściu formalnym ugrzązłem okrutnie.

Celem dzisiejszego wpisu jest pokazanie jak bardzo wektoryzacja może przyspieszyć tego typu symulację. Zaczniemy od podejścia naiwnego:

import numpy as np
import time
czas_startu = time.time()

def symuluj_podziały(ile_poziomów):
    suma_razem = 0
    bieżące_segmenty = [1.0]

    for _ in range(ile_poziomów):
        następne_segmenty = []
        for segment in bieżące_segmenty:
            x = np.random.uniform(0, segment)
            iloczyn = x * (segment - x)
            suma_razem += iloczyn
            następne_segmenty.extend([x, segment - x])
        
        bieżące_segmenty = następne_segmenty
    return suma_razem

wynik = [symuluj_podziały(20) for _ in range(10)]
wynik_średni = np.mean(wynik)

print(wynik_średni)
czas_zakończenia = time.time()
czas_wykonania = czas_zakończenia - czas_startu
print(f"Czas wykonania: {czas_wykonania} sekund")

Na wyjściu dostałem:

0.49985632477171577
Czas wykonania: 24.606609582901 sekund

Jeżeli przyjrzymy się powyższemu kodowi, zobaczymy, że wykonuje on 10 razy funkcję symuluj_podziały(20), która z kolei buduje kolejne podziały odcinków począwszy od pojedynczego odcinka o długości 1 aż do 220=1048576 (milion z niewielkim hakiem), czyli summa summarum musi "przerobić" około dwóch milionów odcinków.

A teraz ten sam pomysł zrealizowany w wersji z wektoryzacją:

import numpy as np
import time
czas_startu = time.time()

def symuluj_podziały(ile_poziomów):
    suma_razem = 0
    bieżące_segmenty = np.array([1.0])

    for poziom in range(ile_poziomów):
        liczba_segmentów = 2 ** (poziom + 1)
        
        x = np.random.uniform(0, 1, liczba_segmentów // 2) * bieżące_segmenty
        
        iloczyny = x * (bieżące_segmenty - x)
        suma_razem += np.sum(iloczyny)
        
        bieżące_segmenty = np.concatenate((x, bieżące_segmenty - x))

    return suma_razem

wyniki = [symuluj_podziały(20) for _ in range(10)]
wynik_średni = np.mean(wyniki)

print(wynik_średni)
czas_zakończenia = time.time()
czas_wykonania = czas_zakończenia - czas_startu
print(f"Czas wykonania: {czas_wykonania} sekund")

Wynik:

0.49983291120363094
Czas wykonania: 0.32900047302246094 sekund

Jak widać identyczny algorytm dzięki zastosowaniu wektoryzacji przyspieszył mniej więcej siedemdziesięciopięciokrotnie! Oznacza to również, że jeżeli dysponujemy budżetem czasowym z oryginalnej wersji skryptu, możemy poprawić dokładność symulacji poprzez zwiększenie liczby rekurencyjnych wywołań z 20 do 26 (bo 26=64):

wyniki = [symuluj_podziały(26) for _ in range(10)]

Wynik:

0.49998513006056455
Czas wykonania: 18.292614221572876 sekund

Możnaby pewnie jeszcze spróbować z PyPy, ale już mi się nie chciało.

https://xpil.eu/gNVza

2 komentarze

  1. Coś oszukujesz. Tu wydruki z mojego programu (czas w sekundach):

    —-
    liczba poziomów: 20
    suma = 0,499853456900863 (1048575 segmentów)
    czas = 00:00:00

    liczba poziomów: 21
    suma = 0,499872360965172 (2097151 segmentów)
    czas = 00:00:00

    liczba poziomów: 22
    suma = 0,499897206843164 (4194303 segmentów)
    czas = 00:00:00

    liczba poziomów: 23
    suma = 0,499918290331034 (8388607 segmentów)
    czas = 00:00:00

    liczba poziomów: 24
    suma = 0,49997957064264 (16777215 segmentów)
    czas = 00:00:01

    liczba poziomów: 25
    suma = 0,499987464149637 (33554431 segmentów)
    czas = 00:00:02

    liczba poziomów: 26
    suma = 0,499994476817779 (67108863 segmentów)
    czas = 00:00:05

    liczba poziomów: 27
    suma = 0,499994351481675 (134217727 segmentów)
    czas = 00:00:11
    —-

    A tu program:

    Procedure symuluj_podzialy(ile_poziomow: integer);
    var i, j, ile_segmentow_poziomu, licznik_segmentow: integer;
        segmenty: array of Extended;
        suma, segmentOld: Extended;
        time1, time2: TDateTime;
    begin
      time1:=GetTime;
      suma:=0;
      licznik_segmentow:=0;
      SetLength(segmenty, 1);
      segmenty[High(segmenty)]:=1.0;
      for i:=1 to ile_poziomow do begin
        ile_segmentow_poziomu:=High(segmenty);
        for j:=0 to ile_segmentow_poziomu do begin
          segmentOld:=segmenty[j];
          segmenty[j]:=Random*segmentOld;
          SetLength(segmenty, Length(segmenty)+1);
          segmenty[High(segmenty)]:=segmentOld-segmenty[j];
          suma:=suma+segmenty[j]*segmenty[High(segmenty)];
          Inc(licznik_segmentow);
        end;
      end;
      time2:=GetTime;
      Form1.Memo1.Lines.Add('liczba poziomów: '+IntToStr(ile_poziomow));
      Form1.Memo1.Lines.Add('suma = '+FloatToStr(suma)+'  ('+IntToStr(licznik_segmentow)+' segmentów)');
      Form1.Memo1.Lines.Add('czas = '+TimeToStr(time2-time1));
    end;
    
    
    procedure TForm1.BStartClick(Sender: TObject);
    begin
      Randomize;
      symuluj_podzialy(27);
    end;
    
    1. Wszystko pięknie, tylko twoja procedura wykonuje obliczenia raz. A moja powtarza symulację dziesięciokrotnie i zwraca średnią:

      ...
      wyniki = [symuluj_podziały(20) for _ in range(10)]
      wynik_średni = np.mean(wyniki)
      ...
      

      Pojedyncze wykonanie z wektoryzacją wykonuje się u mnie w 0.032 sekund, a bez niej w 2.61 sekund (dla 20 poziomów rekurencji).

Leave a Comment

Komentarze mile widziane.

Jeżeli chcesz do komentarza wstawić kod, użyj składni:
[code]
tutaj wstaw swój kod
[/code]

Jeżeli zrobisz literówkę lub zmienisz zdanie, możesz edytować komentarz po jego zatwierdzeniu.