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 2²⁰=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 2⁶=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