V ist (n, p) numpy Array, typischerweise sind die Abmessungen n ~ 10, p ~ 20000

Der Code, den ich jetzt habe, sieht so aus

A = np.zeros(p)
for i in xrange(n):
    for j in xrange(i+1):
        A += F[i,j] * V[i,:] * V[j,:]

Wie würde ich vorgehen, um die doppelte Python-for-Schleife zu vermeiden?

4
user1984528 21 Nov. 2013 im 02:52

3 Antworten

Beste Antwort

Während Isaacs Antwort vielversprechend erscheint, da sie diese beiden verschachtelten for-Schleifen entfernt, müssen Sie ein Zwischenarray M erstellen, das n mal so groß ist wie Ihr ursprüngliches V Array. Python for Loops sind nicht billig, aber der Speicherzugriff ist auch nicht kostenlos:

n = 10
p = 20000
V = np.random.rand(n, p)
F = np.random.rand(n, n)

def op_code(V, F):
    n, p = V.shape
    A = np.zeros(p)
    for i in xrange(n):
        for j in xrange(i+1):
            A += F[i,j] * V[i,:] * V[j,:]
    return A

def isaac_code(V, F):
    n, p = V.shape
    F = F.copy()
    F[np.triu_indices(n, 1)] = 0
    M = (V.reshape(n, 1, p) * V.reshape(1, n, p)) * F.reshape(n, n, 1)
    return M.sum((0, 1))

Wenn Sie jetzt beide für eine Probefahrt nehmen:

In [20]: np.allclose(isaac_code(V, F), op_code(V, F))
Out[20]: True

In [21]: %timeit op_code(V, F)
100 loops, best of 3: 3.18 ms per loop

In [22]: %timeit isaac_code(V, F)
10 loops, best of 3: 24.3 ms per loop

Das Entfernen der for-Schleifen kostet Sie also eine 8-fache Verlangsamung . Keine sehr gute Sache ... An dieser Stelle möchten Sie vielleicht sogar überlegen, ob eine Funktion, deren Auswertung etwa 3 ms dauert, einer weiteren Optimierung bedarf. Falls Sie dies tun, gibt es eine kleine Verbesserung, die mit np.einsum erzielt werden kann:

def einsum_code(V, F):
    n, p = V.shape
    F = F.copy()
    F[np.triu_indices(n, 1)] = 0
    return np.einsum('ij,ik,jk->k', F, V, V)

Und nun:

In [23]: np.allclose(einsum_code(V, F), op_code(V, F))
Out[23]: True

In [24]: %timeit einsum_code(V, F)
100 loops, best of 3: 2.53 ms per loop

Das ist also eine Beschleunigung von ungefähr 20%, die Code einführt, der möglicherweise nicht so lesbar ist wie Ihre for-Schleifen. Ich würde sagen, es lohnt sich nicht ...

10
Jaime 21 Nov. 2013 im 01:35

Das Schwierige daran ist, dass Sie nur die Summe der Elemente mit j <= i nehmen möchten. Wenn nicht, könnten Sie Folgendes tun:

M = (V.reshape(n, 1, p) * V.reshape(1, n, p)) * F.reshape(n, n, 1)
A = M.sum(0).sum(0)

Wenn F symmetrisch ist (wenn F[i,j] == F[j,i]), können Sie die Symmetrie von M wie folgt ausnutzen:

D = M[range(n), range(n)].sum(0)
A = (M.sum(0).sum(0) - D) / 2.0 + D

Das heißt, dies ist wirklich kein großartiger Kandidat für die Vektorisierung, da Sie n << p haben und Ihre for - Schleifen daher keinen großen Einfluss auf die Geschwindigkeit dieser Berechnung haben werden.

Bearbeiten : Wie Bill unten sagte, können Sie nur sicherstellen, dass die Elemente von F, die Sie nicht verwenden möchten, zuerst auf Null und dann auf M.sum(0).sum(0) gesetzt werden. Ergebnis wird sein, was Sie wollen.

7
Isaac 20 Nov. 2013 im 23:22

Der Ausdruck kann geschrieben werden als

formula

Und so können Sie es so mit dem np.newaxis - Konstrukt summieren:

na = np.newaxis
X = (np.tri(n)*F)[:,:,na]*V[:,na,:]*V[na,:,:]
X.sum(axis=1).sum(axis=0)

Hier wird ein 3D-Array X[i,j,p] erstellt, und dann werden die beiden ersten Achsen summiert, was zu einem 1D-Array A[p] führt. Zusätzlich wurde F mit einer dreieckigen Matrix multipliziert, um die Summierung entsprechend dem Problem einzuschränken.

1
flonk 21 Nov. 2013 im 13:22