巨大的对称矩阵——如何巧妙地存储和使用它——Python答案

【问题标题】：Huge symmetric matrix - how to store and use it cleverly - Python巨大的对称矩阵——如何巧妙地存储和使用它——Python
【发布时间】：2015-11-03 15:00:14
【问题描述】：

我有一个对称矩阵。现在，问题是我需要填充这样一个维度矩阵(32**3) x (32**3)。我需要填充矩阵的原因是因为在我的程序中我将它用于各种计算：我正在反转它，我正在将它与其他矩阵相乘......在我看来，为了执行这些不同的您确实需要实际存储完整矩阵的计算，例如，您不能只使用其中的一半（但我可能错了，在这种情况下请告诉我应该怎么做）。

问题是这样的矩阵对于我的计算机来说太大了，我收到以下错误：

Traceback (most recent call last):
    File "program.py", line 191, in <module>
        A = zeros((n_x*n_y*n_z. n_x*n_y*n_z), float)
MemoryError

这里，n_x = 32。那么，我该如何解决这个问题呢？有没有办法存储这么大的矩阵，或者避免存储它的巧妙方法？两种方法对我来说都可以，只要我可以使用它们而不会在计算中出错。

为了完整起见，我在下面报告A 矩阵是如何构建的：

    n_x = n_y = n_z = 32
    L_x = L_y = L_z = n_x
    A = zeros((n_x*n_y*n_z , n_x*n_y*n_z), float)
    P_0 = 50.0
    sigma_x = sigma_y = sigma_z = 0.9
    sigma_root = np.sqrt(sigma_x**2 + sigma_y**2 + sigma_z**2)
    twosigmasquared = 2.*sigma_root**2
    for l in range(n_x*n_y*n_z):
        for m in range(n_x*n_y*n_z):
            A[l][m] = P_0*(L_x/(np.sqrt(2.*np.pi)*sigma_root*n_x**2)) * (L_y/(np.sqrt(2.*np.pi)*sigma_root*n_y**2)) * (L_z/(np.sqrt(2.*np.pi)*sigma_root*n_z**2))*np.exp((-((x[l]-x[m])**2)-((y[l]-y[m])**2)-((z[l]-z[m])**2))/twosigmasquared)
            A[m][l] = A[l][m]

【问题讨论】：

一个明显的解决方案是买一台更大的机器。该矩阵具有~ 10^9 元素，并且应该很容易适应具有~ 16Gb RAM 的机器上的内存。此类机器可在 AWS EC2 上使用。另一种解决方案是动态计算元素：def A(i,j): return ....
取决于numpy 存储矩阵的方式，这可能会消耗大量内存。在它最紧凑的形式（使用对称性）中，它将消耗大约 4Gb 内存，否则它将是 8Gb 内存甚至 16Gb 内存。如果您无法同时在内存中保存整个矩阵，则需要处理该问题。
另一个好问题，别忘了接受对你有更多帮助的答案！

标签： python bigdata matrix-multiplication matrix-inverse symmetric

【解决方案1】：

要节省 50% 的空间，您可以使用 scipy 中的 lil_sparse 矩阵。

from scipy import sparse as S
A = S.lil_matrix((n_x*n_y*n_z , n_x*n_y*n_z), float)

n_x = n_y = n_z = 32
L_x = L_y = L_z = n_x
P_0 = 50.0
sigma_x = sigma_y = sigma_z = 0.9
sigma_root = np.sqrt(sigma_x**2 + sigma_y**2 + sigma_z**2)
twosigmasquared = 2.*sigma_root**2
for l in range(n_x*n_y*n_z):
    for m in range(l, n_x*n_y*n_z): # Filling only the top half
        A[l][m] = P_0*(L_x/(np.sqrt(2.*np.pi)*sigma_root*n_x**2)) * (L_y/(np.sqrt(2.*np.pi)*sigma_root*n_y**2)) * (L_z/(np.sqrt(2.*np.pi)*sigma_root*n_z**2))*np.exp((-((x[l]-x[m])**2)-((y[l]-y[m])**2)-((z[l]-z[m])**2))/twosigmasquared)

那么你可以编写一个辅助函数，而不是访问矩阵本身：

def getA(i, j):
    if i < j:
        return A[j, i]
    else:
        return A[i, j]

但是，如果您想使用标准方法计算矩阵的逆矩阵、想要有效地乘以矩阵或进行任何操作，这对您没有多大帮助。

将整个矩阵保存在内存中可能是更好的选择。

【讨论】：