节省内存的稀疏对称矩阵计算答案

【问题标题】：Memory-efficient sparse symmetric matrix calculations节省内存的稀疏对称矩阵计算
【发布时间】：2018-06-22 19:42:53
【问题描述】：

我必须执行大量这样的计算：

X.dot(Y).dot(Xt)

X = 1 x n 矩阵

Y = 对称 n x n 矩阵，每个元素为 5 个值之一（例如 0、0.25、0.5、0.75、1）

Xt = n x 1 矩阵，X 的转置，即X[np.newaxis].T

X 和 Y 是稠密的。我遇到的问题是大 n，由于内存问题，我无法存储和使用矩阵 Y。我仅限于使用一台机器，所以分布式计算不是一种选择。

我突然想到 Y 有 2 个特性，理论上可以减少存储 Y 所需的内存量：

Y 的元素被一小部分值覆盖。
Y 是对称的。

如何在实践中实现这一点？我查找了对称矩阵的存储，但据我所知，所有 numpy 矩阵乘法都需要“解包”对称性以生成常规的 n x n 矩阵。

我了解 numpy 专为内存计算而设计，因此我愿意接受不依赖于 numpy 的基于 python 的替代解决方案。我也愿意牺牲速度来提高内存效率。

【问题讨论】：

sparse 的用法听起来不适合您的任务（包括标题和标签用法）！此外，它相当广泛。这听起来有点像：我怎样才能为我的任务实现我自己的矩阵例程，这不是微不足道的。特征 1 可以通过使用一个字节或更少作为类型来利用，而特征 2 在稀疏矩阵工具中很常见，所以有很多文献。
@dkato 他说 X 和 Y 是密集的。此外，我认为 scipy 的稀疏矩阵中没有对称性利用。
小 n 的示例 Y：[[0, 0.5, 0.25, 1], [0.5, 0, 0.75, 0.25], [0.25, 0.75, 0, 0.25], [1, 0.25, 0.25, 0]] 我可能对稀疏有误，但我的理解是存储整数索引和对 5 个不同值的引用比为每个元素存储一个浮点数要便宜。我有什么问题吗？
@sascha 对，对不起。我对标题和文字感到困惑。
当然可以更便宜。但是稀疏矩阵通常被定义为（wiki）：在数值分析和计算机科学中，稀疏矩阵或稀疏数组是其中大部分元素为零的矩阵，这是非常不同的。这与存储条目所需的位数无关，而与条目模式的存储方式有关（涉及利用对称性；但与条目的离散性无关）。

标签： python numpy matrix sparse-matrix

【解决方案1】：

更新：我发现使用将 3 个矩阵元素填充到一个字节中的格式实际上非常快。在下面的示例中，与使用@ 的直接乘法相比，速度损失小于2x，而空间节省大于20x。

>>> Y = np.random.randint(0, 5, (3000, 3000), dtype = np.int8)
>>> i, j = np.triu_indices(3000, 1)
>>> Y[i, j] = Y[j, i]
>>> values = np.array([0.3, 0.5, 0.6, 0.9, 2.0])
>>> Ycmp = (np.reshape(Y, (1000, 3, 3000)) * np.array([25, 5, 1], dtype=np.int8)[None, :, None]).sum(axis=1, dtype=np.int8)
>>> 
>>> full = values[Y]
>>> x @ full @ x
1972379.8153972814
>>> 
>>> vtable = values[np.transpose(np.unravel_index(np.arange(125), (5,5,5)))]
>>> np.dot(np.concatenate([(vtable * np.bincount(row, x, minlength=125)[:, None]).sum(axis=0) for row in Ycmp]), x)
1972379.8153972814
>>> 
>>> timeit('x @ full @ x', globals=globals(), number=100)
0.7130507210385986
>>> timeit('np.dot(np.concatenate([(vtable * np.bincount(row, x, minlength=125)[:, None]).sum(axis=0) for row in Ycmp]), x)', globals=globals(), number=100)
1.3755558349657804

以下解决方案速度较慢且内存效率较低。我将它们留作参考。

如果你能负担得起每个矩阵元素半个字节，那么你可以像这样使用np.bincount：

>>> Y = np.random.randint(0, 5, (1000, 1000), dtype = np.int8)
>>> i, j = np.triu_indices(1000, 1)
>>> Y[i, j] = Y[j, i]
>>> values = np.array([0.3, 0.5, 0.6, 0.9, 2.0])
>>> full = values[Y]
>>> # full would correspond to your original matrix,
>>> # Y is the 'compressed' version
>>>
>>> x = np.random.random((1000,))
>>>
>>> # direct method for reference 
>>> x @ full @ x
217515.13954751115
>>> 
>>> # memory saving version
>>> np.dot([(values * np.bincount(row, x)).sum() for row in Y], x)
217515.13954751118
>>>
>>> # to save another almost 50% exploit symmetry
>>> upper = Y[i, j]
>>> diag = np.diagonal(Y)
>>> 
>>> boundaries = np.r_[0, np.cumsum(np.arange(999, 0, -1))]
>>> (values*np.bincount(diag, x*x)).sum() + 2 * np.dot([(values*np.bincount(upper[boundaries[i]:boundaries[i+1]], x[i+1:],minlength=5)).sum() for i in range(999)], x[:-1])
217515.13954751115

【讨论】：

【解决方案2】：

Y 的每一行，如果按照@PaulPanzer 的回答中的建议表示为数据类型为int 的numpy.array，则可以压缩以占用更少的内存：实际上，您可以使用 64 位存储 27 个元素，因为64 / log2(5) = 27.56...

首先，压缩：

import numpy as np

row = np.random.randint(5, size=100)

# pad with zeros to length that is multiple of 27
if len(row)%27:
    row_pad = np.append(row, np.zeros(27 - len(row)%27, dtype=int))
else:
    row_pad = row

row_compr = []
y_compr = 0
for i, y in enumerate(row_pad):
    if i > 0 and i % 27 == 0:
        row_compr.append(y_compr)
        y_compr = 0
    y_compr *= 5
    y_compr += y

# append last 
row_compr.append(y_compr)
row_compr = np.array(row_compr, dtype=np.int64)

然后，解压：

row_decompr = []
for y_compr in row_compr:
    y_block = np.zeros(shape=27, dtype=np.uint8)
    for i in range(27):
        y_block[26-i] = y_compr % 5
        y_compr = int(y_compr // 5)
    row_decompr.append(y_block)

row_decompr = np.array(row_decompr).flatten()[:len(row)]

解压后的数组与Y的原始行重合：

assert np.allclose(row, row_decompr)

【讨论】：