iteración eficiente sobre matriz 3D?

Question

Estoy usando Python y Numpy para hacer algunos análisis de datos.

Tengo una gran matriz 3D (NxNxN), donde cada celda es nuevamente una matriz, esta vez una matriz de 3x3. Llamando a la matriz data, que se ve así:

data[N,N,N,3,3] 

Tengo que encontrar los valores propios de todas las matrices de 3x3, y por eso el uso rutinario de Numpy eigvals, pero lleva mucho tiempo hacer. En este momento, prácticamente hago esto:

for i in range(N): 
    for j in range(N): 
     for k in range(N): 
      a = np.linalg.eigvals(data[i,j,k,:,:]) 

Para N = 256, esto lleva alrededor de una hora. ¿Alguna idea sobre cómo hacer esto más eficiente?

Muchas gracias por cualquier sugerencia!

Answer 1

itertools.product es mejor que los bucles anidados, estéticamente hablando. Pero no creo que haga tu código mucho más rápido. Mi prueba sugiere que la iteración no es su cuello de botella.

>>> bigdata = numpy.arange(256 * 256 * 256 * 3 * 3).reshape(256, 256, 256, 3, 3) 
>>> %timeit numpy.linalg.eigvals(bigdata[100, 100, 100, :, :]) 
10000 loops, best of 3: 52.6 us per loop 

Así subestimar:

>>> .000052 * 256 * 256 * 256/60 
14.540253866666665 

Eso es 14 minutos como mínimo en mi equipo, que es bastante nueva. Veamos cuánto duran los bucles ...

>>> def just_loops(N): 
...  for i in xrange(N): 
...   for j in xrange(N): 
...    for k in xrange(N): 
...     pass 
... 
>>> %timeit just_loops(256) 
1 loops, best of 3: 350 ms per loop 

Órdenes de magnitud menores, como dijo DSM. Incluso el trabajo de cortar la matriz sola es más importante:

>>> def slice_loops(N, data): 
...  for i in xrange(N): 
...   for j in xrange(N): 
...    for k in xrange(N): 
...     data[i, j, k, :, :] 
... 
>>> %timeit slice_loops(256, bigdata) 
1 loops, best of 3: 33.5 s per loop 

Answer 2

Estoy seguro de que hay una buena manera de hacerlo en NumPy, pero en general, itertools.product es más rápido que los bucles anidados en rangos.

from itertools import product 

for i, j, k in product(xrange(N), xrange(N), xrange(N)): 
    a = np.linalg.eigvals(data[i,j,k,:,:]) 

Answer 3

ya que todos los cálculos son independientes, se puede utilizar el módulo de multiprocesamiento para acelerar el cálculo si usted tiene un procesador de varios núcleos.