如何计算一个非常大的相关矩阵

Question

我有 np.array 个观测值 z，其中 z.shape 是 (100000, 60)。我想有效地计算 100000x100000 相关矩阵，然后将这些元素的坐标和值写入磁盘 > 0.95（这只是总数的一小部分）。

我的暴力破解版本如下所示，但不出所料，非常慢：

for i1 in range(z.shape[0]):
    for i2 in range(i1+1):
        r = np.corrcoef(z[i1,:],z[i2,:])[0,1]
        if r > 0.95:
            file.write("%6d %6d %.3f\n" % (i1,i2,r))

我意识到使用 np.corrcoef(z) 可以在一次操作中更有效地计算相关矩阵本身，但是内存需求会很大。我也知道可以将数据集分解成块并一次计算相关矩阵的一小部分，但是对其进行编程和跟踪索引似乎不必要地复杂。

是否有另一种既易于编码又不会对物理内存提出过多要求的方法（例如，使用 memmap 或 pytables）？

Answer 1

根据我的粗略计算，您需要一个包含 100,000^2 个元素的相关矩阵。假设浮动，这将占用大约 40 GB 的内存。
这可能不适合计算机内存，否则您可以使用 corrcoef。有一种基于特征向量的奇特方法，我现在找不到，它进入了（必然）复杂的类别...... 相反，依靠这样一个事实，即对于零均值数据，可以使用点积找到协方差。

z0 = z - mean(z, 1)[:, None]
cov = dot(z0, z0.T)
cov /= z.shape[-1]

这可以通过方差归一化转化为相关性

sigma = std(z, 1)
corr = cov
corr /= sigma
corr /= sigma[:, None]

当然内存占用还是个问题。您可以使用内存映射数组（确保它已打开以进行读写）和 dot 的 out 参数来解决此问题（另一个示例请参见）

N = z.shape[0]
arr = np.memmap('corr_memmap.dat', dtype='float32', mode='w+', shape=(N,N)) 
dot(z0, z0.T, out=arr)
arr /= sigma
arr /= sigma[:, None]

然后你可以遍历结果数组，找到相关系数大的索引。（您可以使用 where(arr > 0.95) 直接找到它们，但是比较会创建一个非常大的布尔数组，它可能适合也可能不适合内存）。

Answer 2

您可以使用 scipy.spatial.distance.pdist 和 metric = correlation 来获得所有没有对称项的相关性。不幸的是，这仍然会给您留下大约 5e10 个术语，这些术语可能会溢出您的记忆。

您可以尝试重新制定 KDTree（理论上可以处理 , and therefore correlation distance) to filter for higher correlations, but with 60 dimensions it's unlikely that would give you much speedup. The curse of dimensionality sucks.

你最好的选择可能是使用 scipy.spatial.distance.cdist(..., metric = correlation) 暴力破解数据块，然后只保留每个块中的高相关性。一旦你知道你的内存可以处理多大的块而不会因为你的计算机的内存架构而变慢，它应该比一次处理一个块快得多。

Answer 3

在尝试了其他人提出的memmap解决方案后，我发现虽然它比我原来的方法（在我的Macbook上花了大约4天）更快，但仍然需要很长时间（至少一天） -- 大概是由于对输出文件的逐个元素写入效率低下。考虑到我需要多次运行计算，这是不可接受的。

最后，（对我而言）最好的解决方案是登录 Amazon Web Services EC2 门户，创建一个 120+ GiB 的虚拟机实例（从配备 Anaconda Python 的映像开始） RAM，上传输入数据文件，并完全在核心内存中进行计算（使用矩阵乘法）。大约两分钟就完成了！

供参考，我用的代码基本是这样的：

import numpy as np
import pickle
import h5py

# read nparray, dimensions (102000, 60)

infile = open(r'file.dat', 'rb')
x = pickle.load(infile)
infile.close()     

# z-normalize the data -- first compute means and standard deviations
xave = np.average(x,axis=1)
xstd = np.std(x,axis=1)

# transpose for the sake of broadcasting (doesn't seem to work otherwise!)
ztrans = x.T - xave
ztrans /= xstd

# transpose back
z = ztrans.T

# compute correlation matrix - shape = (102000, 102000)
arr = np.matmul(z, z.T)   
arr /= z.shape[0]

# output to HDF5 file
with h5py.File('correlation_matrix.h5', 'w') as hf:
        hf.create_dataset("correlation",  data=arr)

Answer 4

请查看 deepgraph 包。

https://deepgraph.readthedocs.io/en/latest/tutorials/pairwise_correlations.html

我试过 z.shape = (2500, 60) 和 pearsonr 2500 * 2500。速度极快。

不确定 100000 x 100000 但值得一试。

如何计算一个非常大的相关矩阵

How to calculate a very large correlation matrix

numpy

large-data

correlation