Beschleunigung für Schleife in Faltung für numpy 3D-Array?

Question

Faltung entlang Z-Vektor eines 3d numpy-Arrays, dann andere Operationen auf den Ergebnissen, aber es ist langsam, wie es jetzt implementiert ist. Ist die for-Schleife, was mich hier verlangsamt oder ist die Faltung? Ich habe versucht, zu einem 1d Vektor umzuformen und die Faltung in 1 Durchlauf durchzuführen (wie ich in Matlab tat), ohne die for Schleife, aber es verbessert nicht Leistung. Meine Matlab-Version ist ungefähr 50% schneller als alles, was ich in Python finden kann. Relevanter Codeabschnitt:

convolved=np.zeros((y_lines,x_lines,z_depth)) 
for i in range(0, y_lines): 
    for j in range(0, x_lines): 
     convolved[i,j,:]= fftconvolve(data[i,j,:], Gauss) #80% of time here 
     result[i,j,:]= other_calculations(convolved[i,j,:]) #20% of time here 

Gibt es einen besseren Weg als eine for-Schleife? Von Cython gehört, aber ich habe jetzt begrenzte Erfahrung in Python, würde auf die einfachste Lösung zielen.

Answer 1

Die von Ihnen verwendete fftconvolve-Funktion stammt vermutlich von SciPy. Wenn ja, beachten Sie, dass N-dimensionale Arrays benötigt werden. So ein schneller Weg, um Ihre Faltung zu tun wäre, um das 3D-Kernel zu erzeugen, die in den x und y Dimensionen und macht eine 1d Gaußsche Faltung in z.

Einige Code und Timing-Ergebnisse sind unten nichts zu tun, entspricht. Auf meinem Rechner und mit einigen Spielzeug Daten, führte dies zu einer 10 × Speedup, wie Sie sehen können:

import numpy as np 
from scipy.signal import fftconvolve 
from scipy.ndimage.filters import gaussian_filter 

# use scipy filtering functions designed to apply kernels to isolate a 1d gaussian kernel 
kernel_base = np.ones(shape=(5)) 
kernel_1d = gaussian_filter(kernel_base, sigma=1, mode='constant') 
kernel_1d = kernel_1d/np.sum(kernel_1d) 

# make the 3d kernel that does gaussian convolution in z axis only 
kernel_3d = np.zeros(shape=(1, 1, 5,)) 
kernel_3d[0, 0, :] = kernel_1d 

# generate random data 
data = np.random.random(size=(50, 50, 50)) 

# define a function for loop based convolution for easy timeit invocation 
def convolve_with_loops(data): 
    nx, ny, nz = data.shape 
    convolved=np.zeros((nx, ny, nz)) 
    for i in range(0, nx): 
     for j in range(0, ny): 
      convolved[i,j,:]= fftconvolve(data[i, j, :], kernel_1d, mode='same') 
    return convolved 

# compute the convolution two diff. ways: with loops (first) or as a 3d convolution (2nd) 
convolved = convolve_with_loops(data) 
convolved_2 = fftconvolve(data, kernel_3d, mode='same') 

# raise an error unless the two computations return equivalent results 
assert np.all(np.isclose(convolved, convolved_2)) 

# time the two routes of the computation 
%timeit convolved = convolve_with_loops(data) 
%timeit convolved_2 = fftconvolve(data, kernel_3d, mode='same') 

---

timeit Ergebnisse:

10 loops, best of 3: 198 ms per loop 
100 loops, best of 3: 18.1 ms per loop 

Answer 2

Ich glaube, Sie haben bereits die source code der fftconvolve Funktion. Normalerweise verwendet es für reale Eingaben die Funktionen numpy.fft.rfftn und .irfftn, die N-dimensionale Transformationen berechnen. Da das Ziel mehrere 1-D-Transformationen zu tun ist, können Sie im Grunde fftconvolve wie folgt umschreiben (vereinfacht):

from scipy.signal.signaltools import _next_regular 

def fftconvolve_1d(in1, in2): 

    outlen = in1.shape[-1] + in2.shape[-1] - 1 
    n = _next_regular(outlen) 

    tr1 = np.fft.rfft(in1, n) 
    tr2 = np.fft.rfft(in2, n) 
    out = np.fft.irfft(tr1 * tr2, n) 

    return out[..., :outlen].copy() 

Und um das gewünschte Ergebnis berechnen:

result = fftconvolve_1d(data, Gauss) 

Das funktioniert, weil numpy.fft.rfft und .irfft (beachten Sie das Fehlen von n im Namen) über eine einzelne Achse des Eingabe-Array (die letzte Achse standardmäßig) zu transformieren. Das ist ungefähr 40% schneller als der OP-Code auf meinem System.

---

Weitere Beschleunigung kann durch Verwendung eines anderen FFT-Back-End erreicht werden.

Zum einen scheinen die Funktionen in scipy.fftpack etwas schneller als ihre Numpy-Äquivalente zu sein. Allerdings ist das Ausgabeformat der Scipy-Varianten ziemlich umständlich (siehe docs) und dies macht es schwierig, die Multiplikation durchzuführen.

Ein anderes mögliches Backend ist FFTW durch den PyFFTW-Wrapper. Nachteile sind, dass Transformationen eine langsame "Planungsphase" vorausgehen und dass die Eingaben 16-Byte-ausgerichtet sein müssen, um die beste Leistung zu erzielen. Dies wird in der pyFFTW tutorial ziemlich gut erklärt. Der resultierende Code könnte zum Beispiel sein:

from scipy.signal.signaltools import _next_regular 
import pyfftw 
pyfftw.interfaces.cache.enable() # Cache for the "planning" 
pyfftw.interfaces.cache.set_keepalive_time(1.0) 

def fftwconvolve_1d(in1, in2): 

    outlen = in1.shape[-1] + in2.shape[-1] - 1 
    n = _next_regular(outlen) 

    tr1 = pyfftw.interfaces.numpy_fft.rfft(in1, n) 
    tr2 = pyfftw.interfaces.numpy_fft.rfft(in2, n) 

    sh = np.broadcast(tr1, tr2).shape 
    dt = np.common_type(tr1, tr2) 
    pr = pyfftw.n_byte_align_empty(sh, 16, dt) 
    np.multiply(tr1, tr2, out=pr) 
    out = pyfftw.interfaces.numpy_fft.irfft(pr, n) 

    return out[..., :outlen].copy() 

mit ausgerichteten Ein- und Cache „Planung“ Ich sah in dem OP den Code eine Beschleunigung von fast 3x über.Speicherausrichtung can be easily checked durch Blick auf die Speicheradresse im ctypes.data Attribut eines Numpy-Array.