Autocorrelación de una matriz multidimensional en numpy

Question

Tengo una matriz bidimensional, es decir, una matriz de secuencias que también son matrices. Para cada secuencia me gustaría calcular la autocorrelación, de modo que para una matriz (5,4), obtendría 5 resultados, o una matriz de dimensión (5,7).

Sé que podría simplemente recorrer la primera dimensión, pero eso es lento y mi último recurso. ¿Hay otra manera?

Gracias!

EDIT:

Sobre la base de la respuesta elegida, más el comentario de MTRW, tengo la siguiente función:

def xcorr(x): 
    """FFT based autocorrelation function, which is faster than numpy.correlate""" 
    # x is supposed to be an array of sequences, of shape (totalelements, length) 
    fftx = fft(x, n=(length*2-1), axis=1) 
    ret = ifft(fftx * np.conjugate(fftx), axis=1) 
    ret = fftshift(ret, axes=1) 
    return ret 

Tenga en cuenta que la longitud es una variable global en mi código, así que asegúrese de declarar eso. Tampoco restringí el resultado a los números reales, ya que también debo tener en cuenta los números complejos.

Answer 1

Usando FFT-based autocorrelation:

import numpy 
from numpy.fft import fft, ifft 

data = numpy.arange(5*4).reshape(5, 4) 
print data 
##[[ 0 1 2 3] 
## [ 4 5 6 7] 
## [ 8 9 10 11] 
## [12 13 14 15] 
## [16 17 18 19]] 
dataFT = fft(data, axis=1) 
dataAC = ifft(dataFT * numpy.conjugate(dataFT), axis=1).real 
print dataAC 
##[[ 14.  8.  6.  8.] 
## [ 126. 120. 118. 120.] 
## [ 366. 360. 358. 360.] 
## [ 734. 728. 726. 728.] 
## [ 1230. 1224. 1222. 1224.]] 

Estoy un poco confundido por su estado de cuenta de la dimensión que tiene la respuesta (5, 7), por lo que tal vez hay algo importante que no estoy entendiendo.

EDIT: A sugerencia de MTRW, una versión acolchada que no todo el documento:

import numpy 
from numpy.fft import fft, ifft 

data = numpy.arange(5*4).reshape(5, 4) 
padding = numpy.zeros((5, 3)) 
dataPadded = numpy.concatenate((data, padding), axis=1) 
print dataPadded 
##[[ 0. 1. 2. 3. 0. 0. 0. 0.] 
## [ 4. 5. 6. 7. 0. 0. 0. 0.] 
## [ 8. 9. 10. 11. 0. 0. 0. 0.] 
## [ 12. 13. 14. 15. 0. 0. 0. 0.] 
## [ 16. 17. 18. 19. 0. 0. 0. 0.]] 
dataFT = fft(dataPadded, axis=1) 
dataAC = ifft(dataFT * numpy.conjugate(dataFT), axis=1).real 
print numpy.round(dataAC, 10)[:, :4] 
##[[ 14.  8.  3.  0.  0.  3.  8.] 
## [ 126. 92. 59. 28. 28. 59. 92.] 
## [ 366. 272. 179. 88. 88. 179. 272.] 
## [ 734. 548. 363. 180. 180. 363. 548.] 
## [ 1230. 920. 611. 304. 304. 611. 920.]] 

Tiene que haber una manera más eficiente de hacer esto, sobre todo porque autocorrelación es simétrica y no hacer aprovecha eso.

Answer 2

Para matrices realmente grandes, es importante tener n = 2 ** p, donde p es un número entero. Esto te ahorrará una gran cantidad de tiempo. Por ejemplo:

def xcorr(x): 
    l = 2 ** int(np.log2(length * 2 - 1)) 
    fftx = fft(x, n = l, axis = 1) 
    ret = ifft(fftx * np.conjugate(fftx), axis = 1) 
    ret = fftshift(ret, axes=1) 
    return ret 

Esto podría darle errores de envolvente. Sin embargo, para matrices grandes la autocorrelación debería ser insignificante cerca de los bordes.

Answer 3

Quizás solo sea una preferencia, pero quería seguir la definición. Personalmente, me resulta un poco más fácil seguirlo de esa manera. Esta es mi implementación para una matriz nd arbitraria.

 
from itertools import product 
from numpy import empty, roll
 

def autocorrelate(x): 
    """ 
    Compute the multidimensional autocorrelation of an nd array. 
    input: an nd array of floats 
    output: an nd array of autocorrelations 
    """ 

    # used for transposes 
    t = roll(range(x.ndim), 1) 

    # pairs of indexes 
    # the first is for the autocorrelation array 
    # the second is the shift 
    ii = [list(enumerate(range(1, s - 1))) for s in x.shape] 

    # initialize the resulting autocorrelation array 
    acor = empty(shape=[len(s0) for s0 in ii]) 

    # iterate over all combinations of directional shifts 
    for i in product(*ii): 
     # extract the indexes for 
     # the autocorrelation array 
     # and original array respectively 
     i1, i2 = asarray(i).T 

     x1 = x.copy() 
     x2 = x.copy() 

     for i0 in i2: 
      # clip the unshifted array at the end 
      x1 = x1[:-i0] 
      # and the shifted array at the beginning 
      x2 = x2[i0:] 

      # prepare to do the same for 
      # the next axis 
      x1 = x1.transpose(t) 
      x2 = x2.transpose(t) 

     # normalize shifted and unshifted arrays 
     x1 -= x1.mean() 
     x1 /= x1.std() 
     x2 -= x2.mean() 
     x2 /= x2.std() 

     # compute the autocorrelation directly 
     # from the definition 
     acor[tuple(i1)] = (x1 * x2).mean() 

    return acor