AGM na sterydach

https://xpil.eu/thj

Osiem lat temu pisałem tu o średniej arytmetyczno-geometrycznej, po naszemu AGM (arithmetic-geometric mean). Tu można poczytać: link 1, link 2.

Okazuje się, że przy odrobinie pomyślunku da się rozciągnąć liczenie takiej średniej na trzy liczby.

Przypomnę najpierw, że AGM dwóch liczb x i y liczymy tak, że najpierw liczymy ich średnią arytmetyczną (czyli suma na pół), potem - geometryczną (czyli iloczyn pod pierwiastek), a to co wyjdzie - czyli dwie nowe liczby - poddajemy tej samej operacji, i tak dalej. Po nieskończenie wielu krokach dostaniemy w wyniku AGM. Przy czym rzecz jasna zamiast nieskończenie wielu kroków wystarczy zrobić tyle kroków, żebyśmy byli zadowoleni z aktualnego przybliżenia. Swoją drogą AGM jest dość szybko zbieżne: po 5 krokach dostajemy wynik z błędem rzędu 10^-12, wystarcza do większości obliczeń.

No dobra, ale to dla dwóch liczb. A co z trzema?

Można zrobić tak: trzy liczby x, y, z dzielimy na trzy pary: (x, y), (x, z), (y, z). Następnie z każdej pary liczymy średnią arytmetyczną (dostajemy w wyniku trzy liczby). Te znów dzielimy na trzy pary, z których z kolei liczymy średnie geometryczne. I tak w kółko, aż się zbiegnie do jednej liczby (z zadaną dokładnością).

Jak by to zapisać w Pythonie?

Dość prosto:

sa = lambda x, y: (x + y) * 0.5 # średnia arytmetyczna
sg = lambda x, y: (x * y) ** 0.5 # średnia geometryczna

x, y, z = 10, 100, 1000

for i in range(30):
    x, y, z = sa(x, y), sa(x, z), sa(y, z)
    x, y, z = sg(x, y), sg(x, z), sg(y, z)

print(x, y, z)

Zbieżność tego procesu jest dużo wolniejsza; dokładność 10^-12 dostajemy po około 25-30 krokach. Ale się zbiega.

Jeżeli ktoś by się temu przyjrzał bliżej, okazałoby się, że AGM liczona dla dwóch zmiennych zbiega się wykładniczo, a dla trzech - liniowo. Ale to już może na kiedy indziej.

A jak by to wyglądało dla większej liczby zmiennych?

Hmmm.

Dla czterech zmiennych mamy sześć możliwych par, a więc w pierwszym kroku dostaniemy sześć liczb (a nie cztery). Z tych sześciu liczb możemy zrobić 15 par i tak dalej. Ucieka nam to wykładniczo. Owszem, gdzieś tam się w końcu zbiegnie, ale dla dłuższych list startowych zabraknie nam pamięci zanim to policzymy.

Ale da się inaczej. Zamiast dzielić n zmiennych na pary, można:

(1) Podzielić je na n różnych grup o liczebności (n-1) każda (takich, że w każdej grupie brakuje innej liczby z oryginalnych n). Na przykład dla czterech zmiennych v, x, y, z dzielimy na grupy (v, x, y),(v, x, z),(v, y, z),(x, y, z).

(2) Dla każdej grupy policzyć średnią arytmetyczną.

(3) Wynik (2) to znów n liczb, powtarzamy manewr z punktu pierwszego...

(4) ... tylko teraz liczymy dla każdej grupy średnią geometryczną.

(5) powtarzamy kroki (1) - (4) aż nam się znudzi.

Efekt?

Sprawdźmy dla czterech zmiennych:

sa = lambda x, y, z: (x + y + z) * (1/3)
sg = lambda x, y, z: (x * y * z) ** (1/3)

v, x, y, z = 10, 100, 1000, 100000

for i in range(27):
    v, x, y, z = sa(x, y, z), sa(v, y, z), sa(v, x, z), sa(v, x, y)
    v, x, y, z = sg(x, y, z), sg(v, y, z), sg(v, x, z), sg(v, x, y)

print(v, x, y, z)

Wynik: 13411.839901919342

Nie wiem jak szybko to się zbiega, ale dokładność 10^-12 dostajemy tu po 18 krokach, więc chyba nieco szybciej niż przy trzech zmiennych.

Dla ogólnego przypadku kod wyglądałby mniej więcej tak:

def sa(x): return sum(x)/len(x)


def sg(x):
    r = 1
    for n in x:
        r = r*n
    return r ** (1/len(x))


x = [10, 200, 3000]

while max(x) != min(x):
    nx1 = []
    for j in range(len(x)):
        tx = x.copy()
        tx.pop(j)
        nx1.append(sa(tx))
    nx2 = []
    for j in range(len(x)):
        tx = nx1.copy()
        tx.pop(j)
        nx2.append(sg(tx))
    x = nx2.copy()

print(x[0])

Wynik: 740.9217572745254

Masz jakieś ciekawsze pomysły na liczenie średnich?

Dzisiejszy wpis sponsorował John.

https://xpil.eu/thj