且构网

分享程序员开发的那些事...
且构网 - 分享程序员编程开发的那些事

用numpy计算向量的自相关

更新时间:2021-10-01 20:42:59

NumPy例程用于一维数组.作为最小"的改进,请对归一化步骤使用向量化操作(在前至后一行中使用np.arange):

The NumPy routines are for 1D arrays. As a "minimal" improvement, use a vectorized operation for the normalisation step (use of np.arange in the before-to-last line):

def vector_autocorrelate(t_array):
  n_vectors = len(t_array)
  # correlate each component indipendently
  acorr = np.array([np.correlate(t_array[:,i],t_array[:,i],'full') for i in xrange(3)])[:,n_vectors-1:]
  # sum the correlations for each component
  acorr = np.sum(acorr, axis = 0)
  # divide by the number of values actually measured and return
  acorr /= (n_vectors - np.arange(n_vectors))
  return acorr

对于更大的数组,应该考虑使用傅里叶变换算法进行关联.如果您对此库感兴趣,请查看库 tidynamics 的示例(免责声明:我写了该库,它仅取决于NumPy.

For larger array sizes, you should consider using the Fourier Transform algorithm for correlation. Check out the examples of the library tidynamics if you are interested in that (disclaimer: I wrote the library, it depends only on NumPy).

作为参考,以下是NumPy代码(我为测试编写的代码),您的代码vector_autocorrelate和tidynamics的时间安排.

For reference, here are the timings for a NumPy-code (that I wrote for testing), your code vector_autocorrelate and tidynamics.

size = [2**i for i in range(4, 17, 2)]

np_time = []
ti_time = []
va_time = []
for s in size:
    data = np.random.random(size=(s, 3))
    t0 = time.time()
    correlate = np.array([np.correlate(data[:,i], data[:,i], mode='full') for i in range(data.shape[1])])[:,:s]
    correlate = np.sum(correlate, axis=0)/(s-np.arange(s))
    np_time.append(time.time()-t0)
    t0 = time.time()
    correlate = tidynamics.acf(data)
    ti_time.append(time.time()-t0)
    t0 = time.time()
    correlate = vector_autocorrelate(data)
    va_time.append(time.time()-t0)

您可以看到结果:

print("size", size)
print("np_time", np_time)
print("va_time", va_time)
print("ti_time", ti_time)

大小[16、64、256、1024、4096、16384、65536]

size [16, 64, 256, 1024, 4096, 16384, 65536]

np_time [0.00023794174194335938,0.0002703666687011719,0.0002713203430175781, 0.001544952392578125、0.0278470516204834、0.36094141006469727、6.922360420227051]

np_time [0.00023794174194335938, 0.0002703666687011719, 0.0002713203430175781, 0.001544952392578125, 0.0278470516204834, 0.36094141006469727, 6.922360420227051]

va_time [0.00021696090698242188,0.0001690387725830078,0.000339508056640625,0.0014629364013671875,0.024930953979492188,0.34442687034606934,7.005630731582642]

va_time [0.00021696090698242188, 0.0001690387725830078, 0.000339508056640625, 0.0014629364013671875, 0.024930953979492188, 0.34442687034606934, 7.005630731582642]

ti_time [0.0011148452758789062,0.0008449554443359375,0.0007512569427490234, 0.0010488033294677734、0.0026645660400390625、0.007939338684082031、0.048232316970825195]

ti_time [0.0011148452758789062, 0.0008449554443359375, 0.0007512569427490234, 0.0010488033294677734, 0.0026645660400390625, 0.007939338684082031, 0.048232316970825195]

或绘制它们

plt.plot(size, np_time)
plt.plot(size, va_time)
plt.plot(size, ti_time)
plt.loglog()

对于非常小的数据序列,"N ** 2"算法不可用.

For anything but very small data series, the "N**2" algorithm is unusable.