具有 numpy 和性能的大型 MATVEC 广播

Question

我正在尝试实现 float64 数据的这种一次性数据转换：

data_new[i,j,k,l,m,:,c] = matvec( data_orig[i,j,k,l,m,:,c], interp_matrix[:,:] )

哪里

"Shape data_orig:     (24,30,12,12,12,1024,2)   [20.4 GB]"
"Shape data_new:      (24,30,12,12,12,7776,2)   [154.8GB]"
"Shape interp_matrix: (1024,7776)"

我的第一个方法已经运行了大约一天，仍然在 c=0 项上。

data_new = np.zeros((24,30,12,12,12,7776,2))
for c in range(num_components):
    print("starting c=%d"%(c))
    data_new[...,c] = numpy.dot( data_orig[...,c],  interp_matrix )

撇开内存问题...（主要是。）我现在可以转移到一个有足够 RAM 的盒子。我需要为未来做一些更强大的事情。

如何在最高 "c" 维度上广播操作？

有了这种能力，我想我可以尝试循环遍历主要维度，并希望启用线程。 numpy.show_config() 告诉我这个 Anaconda 发行版是用 mkl_intel_lp64 和 mkl_intel_thread 构建的。我肯定在第一种方法中只使用一个核心。 Google 搜索让我认为我需要设置 OMP_NUM_THREADS 以便将来尝试启用线程。

这里有一个更简洁的例子：

>> np.dot( np.zeros((3,4,32,2))[...,:,0], np.zeros((32,128)) )
#works. :)

>> np.dot( np.zeros((3,4,32,2)), np.zeros((32,128)) )
ValueError: shapes (3,4,32,2) and (32,128) not aligned: 2 (dim 3) != 32 (dim 0)

Answer 1

在爱因斯坦求和符号中，您要的乘积是：

data_new[i,j,k,l,m,a,c] = data_orig[i,j,k,l,m,b,c] * interp_matrix[b,a]

或者在 numpy 中：

data_new = np.einsum('ijklmbc,ba->ijklmac',data_orig, interp_matrix)

一步一步减少这个，每个都相当于最后一个：

data_new = np.einsum('...bc,ba->...ac', data_orig, interp_matrix)
data_new = np.einsum('...bc,ab->...ac', data_orig, interp_matrix.T)
data_new = np.einsum('ab,...bc->...ac', interp_matrix.T, data_orig)
data_new = np.einsum('...ab,...bc->...ac', interp_matrix.T, data_orig)
                    # ^ This is the summation that matmul performs!
data_new = np.matmul(interp_matrix.T, data_orig)
data_new = interp_matrix.T @ data_orig

具有 numpy 和性能的大型 MATVEC 广播

Large MATVEC broadcasting with numpy and performance

python

numpy

blas