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

【问题标题】：How to calculate a very large correlation matrix如何计算一个非常大的相关矩阵
【发布时间】：2019-02-24 22:21:06
【问题描述】：

我有一个 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）既易于编码又不会对物理内存提出过多要求？

【问题讨论】：

请强调您谈的是内存效率
如果有效，会不会“过于复杂”？
如果有一种 Python 的方式来实现这一点，这将是不必要的复杂，它需要的编码工作量大大减少，因此也不太容易出现编码错误。

标签： numpy correlation large-data

【解决方案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 参数（有关另一个示例，请参见Optimizing my large data code with little RAM）来解决此问题

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) 找到它们，但比较会创建一个非常大的布尔数组，它可能适合也可能不适合内存）。

【讨论】：

我认为计算相关系数并不像取点积那么简单。除了减去平均值之外，您还必须在取点积之前对数据进行 z 归一化（除以每个变量的方差），然后将点积除以样本数 N。
@GrantPetty，我在考虑协方差矩阵。您可以通过额外的步骤对其进行标准化以获取相关矩阵，我已经更新了答案以表明这一点。
我对您使用带有 'r+' 选项的 np.load 感到有些困惑。我不应该使用 arr = np.memmap('corr_memmap.dat', dtype='float32', mode='w+', shape=(N,N))
@GrantPetty，是的，这是正确的做法。如果您有现有文件，load 命令可以工作，但您需要memmap 以及数据类型和形状来创建文件。感谢您了解这一点。
我发现相对于我原来的方法，memmap 方法并没有明显帮助加快速度——24 小时后，程序仍在运行，大概是因为它不断地撞击磁盘执行矩阵乘法。

【解决方案2】：

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

您可以尝试重新制定KDTree（理论上可以处理cosine distance，因此可以处理相关距离）来过滤更高的相关性，但是对于 60 个维度，它不太可能给您带来很大的加速。 The curse of dimensionality sucks.

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

【讨论】：

【解决方案3】：

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

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

作为参考，我使用的代码基本上是这样的：

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)

【讨论】：

【解决方案4】：

请查看 deepgraph 包。

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

我在 z.shape = (2500, 60) 和 pearsonr 上尝试了 2500 * 2500。它的速度非常快。

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

【讨论】：