PyFFTW perfomance su array multidimensionali

Question

Ho un array nD, diciamo di dimensioni: (144, 522720) e ho bisogno di calcolare la sua FFT.

PyFFTWsembra più lento di numpye scipy, che NON è previsto.

Sto facendo qualcosa di evidentemente sbagliato?

Di seguito è il mio codice

import numpy
import scipy      
import pyfftw
import time

n1 = 144
n2 = 522720
loops = 2

pyfftw.config.NUM_THREADS = 4
pyfftw.config.PLANNER_EFFORT = 'FFTW_ESTIMATE'
# pyfftw.config.PLANNER_EFFORT = 'FFTW_MEASURE'

Q_1 = pyfftw.empty_aligned([n1, n2], dtype='float64')
Q_2 = pyfftw.empty_aligned([n1, n2], dtype='complex_')
Q_ref = pyfftw.empty_aligned([n1, n2], dtype='complex_')

# repeat a few times to see if pyfft planner helps
for i in range(0,loops):
    Q_1 = numpy.random.rand(n1,n2)

    s1 = time.time()
    Q_ref = numpy.fft.fft(Q_1, axis=0)
    print('NUMPY - elapsed time: ', time.time() - s1, 's.')

    s1 = time.time()
    Q_2 = scipy.fft.fft(Q_1, axis=0)
    print('SCIPY - elapsed time: ', time.time() - s1, 's.')
    print('Equal = ', numpy.allclose(Q_2, Q_ref))

    s1 = time.time()
    Q_2 = pyfftw.interfaces.numpy_fft.fft(Q_1, axis=0)
    print('PYFFTW NUMPY - elapsed time = ', time.time() - s1, 's.')
    print('Equal = ', numpy.allclose(Q_2, Q_ref))

    s1 = time.time()
    Q_2 = pyfftw.interfaces.scipy_fftpack.fft(Q_1, axis=0)
    print('PYFFTW SCIPY - elapsed time = ', time.time() - s1, 's.')
    print('Equal = ', numpy.allclose(Q_2, Q_ref))

    s1 = time.time()
    fft_object = pyfftw.builders.fft(Q_1, axis=0)
    Q_2 = fft_object()
    print('FFTW PURE Elapsed time = ', time.time() - s1, 's')
    print('Equal = ', numpy.allclose(Q_2, Q_ref))

Answer 1

Innanzitutto, se attivi la cache prima del ciclo principale, le interfacce funzionano in gran parte come previsto:

pyfftw.interfaces.cache.enable()
pyfftw.interfaces.cache.set_keepalive_time(30)

È interessante notare che, nonostante la saggezza che dovrebbe essere immagazzinata, la costruzione degli pyfftwoggetti è ancora piuttosto lenta quando la cache è disattivata. Non importa, questo è esattamente lo scopo della cache. Nel tuo caso devi rendere il tempo di mantenimento della cache piuttosto lungo perché il tuo ciclo è molto lungo.

In secondo luogo, non è un confronto equo includere il tempo di costruzione del fft_objectnel test finale. Se lo sposti fuori dal timer, la chiamata fft_objectè una misura migliore.

In terzo luogo, è anche interessante vedere che anche con la cache attivata, la chiamata a numpy_fftè più lenta della chiamata a scipy_fft. Poiché non vi è alcuna differenza evidente nel percorso del codice, suggerisco che si tratti di un problema di memorizzazione nella cache. Questo è il tipo di problema che timeitcerca di mitigare. Ecco il mio codice di temporizzazione proposto che è più significativo:

import numpy
import scipy
import pyfftw
import timeit

n1 = 144
n2 = 522720

pyfftw.config.NUM_THREADS = 4
pyfftw.config.PLANNER_EFFORT = 'FFTW_MEASURE'

Q_1 = pyfftw.empty_aligned([n1, n2], dtype='float64')

pyfftw.interfaces.cache.enable()
pyfftw.interfaces.cache.set_keepalive_time(30)

times = timeit.repeat(lambda: numpy.fft.fft(Q_1, axis=0), repeat=5, number=1)
print('NUMPY fastest time = ', min(times))

times = timeit.repeat(lambda: scipy.fft.fft(Q_1, axis=0), repeat=5, number=1)
print('SCIPY fastest time = ', min(times))

times = timeit.repeat(
    lambda: pyfftw.interfaces.numpy_fft.fft(Q_1, axis=0), repeat=5, number=1)
print('PYFFTW NUMPY fastest time = ', min(times))

times = timeit.repeat(
    lambda: pyfftw.interfaces.scipy_fftpack.fft(Q_1, axis=0), repeat=5, number=1)
print('PYFFTW SCIPY fastest time = ', min(times))

fft_object = pyfftw.builders.fft(Q_1, axis=0)
times = timeit.repeat(lambda: fft_object(Q_1), repeat=5, number=1)
print('FFTW PURE fastest time = ', min(times))

Sulla mia macchina questo dà un output come:

NUMPY fastest time =  0.6622681759763509
SCIPY fastest time =  0.6572431400418282
PYFFTW NUMPY fastest time =  0.4003451430471614
PYFFTW SCIPY fastest time =  0.40362057799939066
FFTW PURE fastest time =  0.324020683998242

Puoi fare un po 'meglio se non lo costringi a copiare l'array di input in un tipo di dati complesso cambiando Q_1in complex128:

NUMPY fastest time =  0.6483533839927986
SCIPY fastest time =  0.847397351055406
PYFFTW NUMPY fastest time =  0.3237176960101351
PYFFTW SCIPY fastest time =  0.3199474769644439
FFTW PURE fastest time =  0.2546963169006631

Questo interessante scipyrallentamento è ripetibile.

Detto questo, se il tuo input è reale, dovresti eseguire una trasformazione reale (per> 50% di velocità con pyfftw) e manipolare l'output complesso risultante.

La cosa interessante di questo esempio è (credo) quanto sia importante la cache nei risultati (che suggerisco è il motivo per cui passare a una trasformazione reale è così efficace nell'accelerare le cose). Vedi qualcosa di drammatico anche quando usi la modifica della dimensione dell'array in 524288 (la potenza successiva di due, che pensi potrebbe forse accelerare le cose, ma non rallentarla drasticamente). In questo caso tutto rallenta un po ', scipysoprattutto. Mi sembra che scipysia più sensibile alla cache, il che spiegherebbe il rallentamento con la modifica dell'input in complex128(522720 è un numero abbastanza carino per FFTing, quindi forse dovremmo aspettarci un rallentamento).

Infine, se la velocità è secondaria alla precisione, puoi sempre utilizzare float a 32 bit come tipo di dati. Se lo combini con una trasformazione reale, ottieni un fattore di accelerazione migliore di 10 rispetto al numpymigliore iniziale dato sopra:

PYFFTW NUMPY fastest time =  0.09026529802940786
PYFFTW SCIPY fastest time =  0.1701313250232488
FFTW PURE fastest time =  0.06202622700948268

(numpy e scipy non cambiano molto perché penso che utilizzino internamente float a 64 bit).

Modifica: ho dimenticato che i fftpackveri FFT di Scipy hanno una strana struttura di output, che si pyfftwreplica con qualche rallentamento. Questo è cambiato per essere più sensato nel nuovo modulo FFT .

La nuova interfaccia FFT è implementata in pyFFTW e dovrebbe essere preferita. Sfortunatamente c'è stato un problema con i documenti in fase di ricostruzione, quindi i documenti erano obsoleti da molto tempo e non mostravano la nuova interfaccia - si spera che ora sia stato risolto.