webentwicklung-frage-antwort-db.com.de

Wenden Sie auf jede Zeile eines Narrays eine Funktion an

Ich habe diese Funktion, um den quadrierten Mahalanobis-Abstand des Vektors x zu berechnen:

def mahalanobis_sqdist(x, mean, Sigma):
   '''
    Calculates squared Mahalanobis Distance of vector x 
    to distibutions' mean 
   '''
   Sigma_inv = np.linalg.inv(Sigma)
   xdiff = x - mean
   sqmdist = np.dot(np.dot(xdiff, Sigma_inv), xdiff)
   return sqmdist

Ich habe ein numpy-Array, das die Form von (25, 4) hat. Ich möchte diese Funktion also auf alle 25 Zeilen meines Arrays ohne for-Schleife anwenden. Wie kann ich also die vektorisierte Form dieser Schleife schreiben?

for r in d1:
    mahalanobis_sqdist(r[0:4], mean1, Sig1)

wo mean1 und Sig1 sind:

>>> mean1
array([ 5.028,  3.48 ,  1.46 ,  0.248])
>>> Sig1 = np.cov(d1[0:25, 0:4].T)
>>> Sig1
array([[ 0.16043333,  0.11808333,  0.02408333,  0.01943333],
       [ 0.11808333,  0.13583333,  0.00625   ,  0.02225   ],
       [ 0.02408333,  0.00625   ,  0.03916667,  0.00658333],
       [ 0.01943333,  0.02225   ,  0.00658333,  0.01093333]])

Ich habe folgendes versucht, aber es hat nicht funktioniert:

>>> vecdist = np.vectorize(mahalanobis_sqdist)
>>> vecdist(d1, mean1, Sig1)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python2.7/dist-packages/numpy/lib/function_base.py", line 1862, in __call__
    theout = self.thefunc(*newargs)
  File "<stdin>", line 6, in mahalanobis_sqdist
  File "/usr/lib/python2.7/dist-packages/numpy/linalg/linalg.py", line 445, in inv
    return wrap(solve(a, identity(a.shape[0], dtype=a.dtype)))
IndexError: Tuple index out of range
13
Vahid Mir

Um eine Funktion auf jede Zeile eines Arrays anzuwenden, können Sie Folgendes verwenden:

np.apply_along_axis(mahalanobis_sqdist, 1, d1, mean1, Sig1)    

In diesem Fall gibt es jedoch einen besseren Weg. Sie müssen nicht auf jede Zeile eine Funktion anwenden. Sie können stattdessen NumPy-Operationen auf das gesamte d1-Array anwenden, um dasselbe Ergebnis zu berechnen. np.einsum kann den for-loop und die beiden Aufrufe von np.dot ersetzen:


def mahalanobis_sqdist2(d, mean, Sigma):
   Sigma_inv = np.linalg.inv(Sigma)
   xdiff = d - mean
   return np.einsum('ij,im,mj->i', xdiff, xdiff, Sigma_inv)

Hier sind einige Benchmarks:

import numpy as np
np.random.seed(1)

def mahalanobis_sqdist(x, mean, Sigma):
   '''
   Calculates squared Mahalanobis Distance of vector x 
   to distibutions mean 
   '''
   Sigma_inv = np.linalg.inv(Sigma)
   xdiff = x - mean
   sqmdist = np.dot(np.dot(xdiff, Sigma_inv), xdiff)
   return sqmdist

def mahalanobis_sqdist2(d, mean, Sigma):
   Sigma_inv = np.linalg.inv(Sigma)
   xdiff = d - mean
   return np.einsum('ij,im,mj->i', xdiff, xdiff, Sigma_inv)

def using_loop(d1, mean, Sigma):
    expected = []
    for r in d1:
        expected.append(mahalanobis_sqdist(r[0:4], mean1, Sig1))
    return np.array(expected)

d1 = np.random.random((25,4))
mean1 = np.array([ 5.028,  3.48 ,  1.46 ,  0.248])
Sig1 = np.cov(d1[0:25, 0:4].T)

expected = using_loop(d1, mean1, Sig1)
result = np.apply_along_axis(mahalanobis_sqdist, 1, d1, mean1, Sig1)
result2 = mahalanobis_sqdist2(d1, mean1, Sig1)
assert np.allclose(expected, result)
assert np.allclose(expected, result2)

In [92]: %timeit mahalanobis_sqdist2(d1, mean1, Sig1)
10000 loops, best of 3: 31.1 µs per loop
In [94]: %timeit using_loop(d1, mean1, Sig1)
1000 loops, best of 3: 569 µs per loop
In [91]: %timeit np.apply_along_axis(mahalanobis_sqdist, 1, d1, mean1, Sig1)
1000 loops, best of 3: 806 µs per loop

Damit ist mahalanobis_sqdist2 etwa 18x schneller als ein for-loop und 26x schneller als die Verwendung von np.apply_along_axis.


Beachten Sie, dass np.apply_along_axis, np.vectorize, np.frompyfunc Python-Utility-Funktionen sind. Unter der Haube verwenden sie for- oder while-loops. Hier gibt es keine echte "Vektorisierung". Sie können syntaktische Hilfestellung leisten, erwarten jedoch nicht, dass Ihr Code besser funktioniert als ein for-loop, den Sie selbst schreiben.

20
unutbu

Die Antwort von @unutbu funktioniert sehr gut, um jede Funktion auf die Zeilen eines Arrays anzuwenden. In diesem speziellen Fall gibt es einige mathematische Symmetrien, die Sie verwenden können, was die Arbeit erheblich beschleunigt, wenn Sie mit großen Arrays arbeiten .

Hier ist eine modifizierte Version Ihrer Funktion:

def mahalanobis_sqdist3(x, mean, Sigma):
    Sigma_inv = np.linalg.inv(Sigma)
    xdiff = x - mean
    return (xdiff.dot(Sigma_inv)*xdiff).sum(axis=-1)

Wenn Sie am Ende irgendeine Art von Sigmaverwenden, würde ich empfehlen, dass Sie Sigma_inv zwischenspeichern und stattdessen als Argument an Ihre Funktion übergeben. Da es in diesem Beispiel 4x4 ist, spielt dies keine Rolle. Ich werde sowieso jedem zeigen, der mit großen Sigmaumgeht.

Wenn Sie nicht dasselbe Sigmanicht wiederholt verwenden, können Sie es nicht zwischenspeichern. Statt die Matrix zu invertieren, können Sie das lineare System mit einer anderen Methode lösen. Hier Ich benutze die in SciPy eingebaute LU -Zerlegung. Dies verbessert die Zeit nur, wenn die Anzahl der Spalten von xrelativ zu der Anzahl der Zeilen groß ist.

Hier ist eine Funktion, die diesen Ansatz zeigt:

from scipy.linalg import lu_factor, lu_solve
def mahalanobis_sqdist4(x, mean, Sigma):
    xdiff = x - mean
    Sigma_inv = lu_factor(Sigma)
    return (xdiff.T*lu_solve(Sigma_inv, xdiff.T)).sum(axis=0)

Hier sind einige Timings. Ich füge die Version mit einsumhinzu, wie in der anderen Antwort erwähnt.

import numpy as np
Sig1 = np.array([[ 0.16043333,  0.11808333,  0.02408333,  0.01943333],
                 [ 0.11808333,  0.13583333,  0.00625   ,  0.02225   ],
                 [ 0.02408333,  0.00625   ,  0.03916667,  0.00658333],
                 [ 0.01943333,  0.02225   ,  0.00658333,  0.01093333]])
mean1 = np.array([ 5.028,  3.48 ,  1.46 ,  0.248])
x = np.random.Rand(25, 4)
%timeit np.apply_along_axis(mahalanobis_sqdist, 1, x, mean1, Sig1)
%timeit mahalanobis_sqdist2(x, mean1, Sig1)
%timeit mahalanobis_sqdist3(x, mean1, Sig1)
%timeit mahalanobis_sqdist4(x, mean1, Sig1)

geben:

1000 loops, best of 3: 973 µs per loop
10000 loops, best of 3: 36.2 µs per loop
10000 loops, best of 3: 40.8 µs per loop
10000 loops, best of 3: 83.2 µs per loop

Wenn Sie jedoch die Größe der betroffenen Arrays ändern, werden die Timing-Ergebnisse geändert. Wenn Sie beispielsweise x = np.random.Rand(2500, 4) angeben, sind die Timings:

10 loops, best of 3: 95 ms per loop
1000 loops, best of 3: 355 µs per loop
10000 loops, best of 3: 131 µs per loop
1000 loops, best of 3: 337 µs per loop

Und wenn wir x = np.random.Rand(1000, 1000), Sigma1 = np.random.Rand(1000, 1000) und mean1 = np.random.Rand(1000) lassen, sind die Timings:

1 loops, best of 3: 1min 24s per loop
1 loops, best of 3: 2.39 s per loop
10 loops, best of 3: 155 ms per loop
10 loops, best of 3: 99.9 ms per loop

Edit : Ich habe bemerkt, dass eine der anderen Antworten die Cholesky-Zerlegung verwendet hat. Da Sigmasymmetrisch und positiv bestimmt ist, können wir tatsächlich bessere Ergebnisse erzielen als meine obigen Ergebnisse. Es gibt Einige gute Routinen von BLAS und LAPACK, die über SciPy erhältlich sind, können mit symmetrisch positiv definierten Matrizen arbeiten. Hier sind zwei schnellere Versionen.

from scipy.linalg.fblas import dsymm
def mahalanobis_sqdist5(x, mean, Sigma_inv):
    xdiff = x - mean
    Sigma_inv = la.inv(Sigma)
    return np.einsum('...i,...i->...',dsymm(1., Sigma_inv, xdiff.T).T, xdiff)
from scipy.linalg.flapack import dposv
def mahalanobis_sqdist6(x, mean, Sigma):
    xdiff = x - mean
    return np.einsum('...i,...i->...', xdiff, dposv(Sigma, xdiff.T)[1].T)

Der erste invertiert immer noch Sigma. Wenn Sie das Inverse vorberechnen und erneut verwenden, ist es wesentlich schneller (der 1000x1000-Fall nimmt auf meinem Computer 35,6 ms mit dem vorberechneten Inverse). Ich habe auch einsum verwendet, um das Produkt dann auf der letzten Achse zusammenzufassen. Dies wurde am Ende etwas schneller als etwas wie (A * B).sum(axis=-1). . Diese beiden Funktionen geben die folgenden Timings ab:

Erster Testfall:

10000 loops, best of 3: 55.3 µs per loop
100000 loops, best of 3: 14.2 µs per loop

Zweiter Testfall:

10000 loops, best of 3: 121 µs per loop
10000 loops, best of 3: 79 µs per loop

Dritter Testfall:

10 loops, best of 3: 92.5 ms per loop
10 loops, best of 3: 48.2 ms per loop
8
IanH

Habe gerade einen wirklich netten Kommentar zu reddit gesehen, der die Dinge vielleicht noch etwas beschleunigen könnte:

Dies ist für niemanden überraschend, der regelmäßig Numpy verwendet. Denn Schleifen In Python sind furchtbar langsam. Eigentlich ist einsum auch ziemlich langsam. Hier ist eine Version, die schneller ist, wenn Sie viele Vektoren haben (500 Vektoren in 4 Dimensionen reichen aus, um diese Version schneller zu machen als Einsum auf meiner Maschine):

def no_einsum(d, mean, Sigma):
    L_inv = np.linalg.inv(numpy.linalg.cholesky(Sigma))
    xdiff = d - mean
    return np.sum(np.dot(xdiff, L_inv.T)**2, axis=1)

Wenn Ihre Punkte auch hochdimensional sind, ist die Berechnung der Umkehrung Langsam (und im Allgemeinen ohnehin eine schlechte Idee), und Sie können Zeit sparen, indem Sie Das System direkt lösen (500 Vektoren in 250 Dimensionen sind ausreichend). .____.] um diese Version auf meinem Rechner am schnellsten zu machen):

def no_einsum_solve(d, mean, Sigma):
    L = numpy.linalg.cholesky(Sigma)
    xdiff = d - mean
    return np.sum(np.linalg.solve(L, xdiff.T)**2, axis=0)
5
Sebastian

Das Problem ist, dass np.vectorize über alle Argumente vektorisiert, aber Sie müssen nur über das erste vektorisieren. Sie müssen das excluded-Schlüsselwortargument für vectorize verwenden:

np.vectorize(mahalanobis_sqdist, excluded=[1, 2])
0
abacabadabacaba