在 python 中为大型 scipy.sparse 矩阵运算分配内存答案

【问题标题】：allocate memory in python for large scipy.sparse matrix operations在 python 中为大型 scipy.sparse 矩阵运算分配内存
【发布时间】：2015-12-22 17:07:05
【问题描述】：

有没有办法为 scipy 稀疏矩阵函数分配内存来处理大型数据集？

具体来说，我正在尝试使用非对称最小二乘平滑（翻译成 python here 和原始 here）对大型质谱数据集（长度约为 60,000）执行基线校正。

该函数（见下文）使用 scipy.sparse 矩阵运算。

def baseline_als(y, lam, p, niter):
  L = len(y)
  D = sparse.csc_matrix(np.diff(np.eye(L), 2))
  w = np.ones(L)
  for i in xrange(niter):
    W = sparse.spdiags(w, 0, L, L)
    Z = W + lam * D.dot(D.transpose())
    z = spsolve(Z, w*y)
    w = p * (y > z) + (1-p) * (y < z)
  return z

当我传递长度为 10,000 或更少的数据集时，我没有问题：

baseline_als(np.ones(10000),100,0.1,10)

但是当传递更大的数据集时，例如

baseline_als(np.ones(50000), 100, 0.1, 10)

我得到一个内存错误，对于该行

  D = sparse.csc_matrix(np.diff(np.eye(L), 2))

【问题讨论】：

不，您不能为 numpy 数组分配内存。我怀疑您在np.eye 或np.diff 函数中遇到了内存错误。两者都将产生(L,L) 形状数组 - 非常大的数组。试试那些没有sparse 电话的人。
@hpaulj 你是对的。当 np.diff 函数传递大 np.eye 生成的数组时出现错误： np.diff(np.eye(len(np.ones(50000)))) 我是 python 新手。有人知道解决方法吗？
通常当人们制作大型稀疏矩阵时，他们会尝试构建它们而不首先制作等效的密集数组。在小情况下使用密集的很方便，但会破坏使用稀疏的许多优点。
我也在 python 中实现了 ALS，我发现使用 solve_banded 可以显着加快速度。如果您愿意，我也愿意分享我的实现。
@perimosocordiae 如果您能分享实现，我将不胜感激。我计划分析一个相当大的数据集，任何速度上的改进都是有益的。

标签： python memory scipy ipython sparse-matrix

【解决方案1】：

尝试改变

D = sparse.csc_matrix(np.diff(np.eye(L), 2))

到

diag = np.ones(L - 2)
D = sparse.spdiags([diag, -2*diag, diag], [0, -1, -2], L, L-2)

D 将是DIAgonal 格式的稀疏矩阵。如果发现 CSC 格式很重要，请使用tocsc() 方法进行转换：

D = sparse.spdiags([diag, -2*diag, diag], [0, -1, -2], L, L-2).tocsc()

以下示例显示新旧版本生成相同的矩阵：

In [67]: from scipy import sparse

In [68]: L = 8

原文：

In [69]: D = sparse.csc_matrix(np.diff(np.eye(L), 2))

In [70]: D.A
Out[70]: 
array([[ 1.,  0.,  0.,  0.,  0.,  0.],
       [-2.,  1.,  0.,  0.,  0.,  0.],
       [ 1., -2.,  1.,  0.,  0.,  0.],
       [ 0.,  1., -2.,  1.,  0.,  0.],
       [ 0.,  0.,  1., -2.,  1.,  0.],
       [ 0.,  0.,  0.,  1., -2.,  1.],
       [ 0.,  0.,  0.,  0.,  1., -2.],
       [ 0.,  0.,  0.,  0.,  0.,  1.]])

新版本：

In [71]: diag = np.ones(L - 2)

In [72]: D = sparse.spdiags([diag, -2*diag, diag], [0, -1, -2], L, L-2)

In [73]: D.A
Out[73]: 
array([[ 1.,  0.,  0.,  0.,  0.,  0.],
       [-2.,  1.,  0.,  0.,  0.,  0.],
       [ 1., -2.,  1.,  0.,  0.,  0.],
       [ 0.,  1., -2.,  1.,  0.,  0.],
       [ 0.,  0.,  1., -2.,  1.,  0.],
       [ 0.,  0.,  0.,  1., -2.,  1.],
       [ 0.,  0.,  0.,  0.,  1., -2.],
       [ 0.,  0.,  0.,  0.,  0.,  1.]])

【讨论】：