Python: Schneller als man denkt - Das Beispiel

In der vorigen Ausgabe habe ich eingeführt, warum Python potenziell langsam ist. Dies in einem Beispiel einmal darzustellen ist der nächste Schritt. Ich werde darüber schreiben, weshalb ich das Beispiel gewählt habe, wie die (reine Python) Implementierung aussieht.

Das Beispiel, welches wir verwenden soll:

• Einem realistischen Szenario entsprechen

• Mehrstufiges Optimierungspotential aufweisen

Für mich fällt die Wahl daher auf die Kernel Density Estimation, ein nicht parametrisches Verfahren zur Schätzung der Dichtefunktion einer Stichprobe. Die Schätzung der Dichtefunktion der Stichprobe nimmt dabei eine Verteilungsfunktion an, welche der Stichprobe folgt. Dieser "Kernel" ist beliebig austauschbar und beeinflusst maßgeblich, wie sich die Dichtefunktion an den Rändern der Stichprobe verhält.

Standardmäßig wird die Dichtefunktion der Gaussverteilung verwendet, sie entsteht durch die Überlagerung der standardverteilte Gaussfunktion jedes Punktes aus der Stichprobe.

$$
p(\mathbf{x}, \mathbf{x'}) =
\frac{1}{N} \sum_{n=1}^{N}
\frac{1}{\sqrt{2\pi h^2}}
\exp\left(
-\frac{\lVert \mathbf{x'} - \mathbf{x}_n \rVert^2}{2h^2}
\right)
$$

Die Formel beschreibt, wie wir diese Dichtefunktion finden:

• $\mathbf{x}$: Die Stichpunkte, welche wir versuchen mit der Dichtefunktion zu beschreiben.

• $\mathbf{x'}$: An neuen Stichtpunkten taste ich die Dichtefunktion ab, welchen Wert bekommen wir für z.B. 3,99?

• $h$: Nimmt den Wert der Standardverteilung an ($h = 1$), je kleiner der Wert, desto stärker folgt die Dichtefunktion Änderungen in der Stichprobe.

• $N$: Anzahl der Punkte in der Stichprobe

• $||\mathbf{x}||$: Die euklidische Norm, berechnet sich aus $\sqrt{x_1^2 + x_2^2 + \dots + x_m^2}$ für einen Vektor der Länge $m$.

Implementierung

Untersuchen wir die Funktion, so gibt es zwei Schleifen:

• $\sum^N_{n=1} \dots$ fordert von uns über alle Elemente der Stichprobe zu iterieren und die darin enthaltenen Term zu berechnen. Da es links in der Funktion steht, ist dies unsere äußere Schleife.

• $|| \mathbf{x'} - \mathbf{x}_n||^2$ fordert uns auf, das wir die quadrierte euklidische Norm für jedes Element des Vektors $\mathbf{x'}$ abzüglich des $n$ten Elements des Vektors $\mathbf{x}$ berechnen. Die Quadrierung und Wurzel sind Invers und können damit gestrichen werden, es bleibt die Summe der einzelnen Quadrate des Distanzmaßes übrig.

Den Term $\sqrt{2\pi h^2}$ ziehen wir aus der Summe raus, da wir ihn nur einmal berechnen müssen. Zur Übersichtlichkeit wird der verbliebene quadratische Term aus der euklidischen Norm als eigene Variable gespeichert.

Die Funktion berechnet für jedes abgetastete $x$ aus allen Datenpunkten $data$ der Stichprobe die Dichte.

def kde_gauss_v1(x: list, data: list, h: float = 1.0):  
    res = []
    term_1 = 1.0/sqrt(2.0 * pi * h**2)

    for _x in x:
        res.append(0.0)
        for d in data:
            term_2 = - (_x - d)**2/(2.0*h**2)
            res[-1] += term_1*exp(term_2)
        res[-1] = 1/len(data)*res[-1]
    return res

Berechnen wir mit dieser Funktion die Dichtefunktion für $M$ Datenpunkte und einer Stichprobe mit $N$ Datenpunkte, so ist die Laufzeit polynomial $O(M \cdot N)$.

Für die Generierung des Titelbildes aus vier verschiedenen Normalverteilungen mit jeweils 2.000 Datenpunkten und einem äquidistanten Raster aus 10.000 abgetasteten Werten der Dichtefunktion brauchte der Rechner 23,2 Sekunden.

Die Laufzeit wurde über 100 Messungen ermittelt, für Python typisch benötigte die erste Runde länger, da hier der Bytecode der Funktion generiert werden musste.