Scipy稀疏特征求解器：循环多次后内存错误且未写入新内容

Question

我正在使用Python和Scipy来对稀疏矩阵进行对角化，这些矩阵的对角线元素是随机的；特别是，我需要获取谱中间的特征值。我写的代码已经运行了几个月，但现在我在处理更大的矩阵时遇到了“内存错误”。让我感到困惑和抓狂的是，这个错误只在构建随机矩阵和对角化的第9次迭代后出现，但我看不出我的代码在每次迭代中存储了额外的内存，因此无法理解为什么在第9次迭代时会失败，而在第1次时却没有。

以下是一些细节（如果我遗漏了什么，提前道歉，我还是新手）：

我构建的每个矩阵都是16000x16000的，里面有15x16000个非零元素。当我处理4000x4000大小的矩阵时，一切都运行得很好。我的代码大部分是

#Initialization
#...

for i in range(dim):
    for n in range(N):
        digit = (i % 2**(n+1)) / 2**n
        index = (i % 2**n) + ((digit + 1) % 2)*(2**n) + (i / 2**(n+1))*(2**(n+1))
        row[dim + N*i + n] = index
        col[dim + N*i + n] = i
        dat[dim + N*i + n] = -G

e_list = open(e_list_name + "_%03dk_%010ds" % (num_states, int(start_time)), "w")
e_log = open(e_log_name + "_%03dk_%010ds" % (num_states, int(start_time)), "w")

for t in range(num_itr): #Begin iterations
    dat[0:dim] = math.sqrt(N/2.0)*np.random.randn(dim) #Get new diagonal elements
    H = sparse.csr_matrix((dat, (row, col))) #Construct new matrix
    vals = sparse.linalg.eigsh(H, k = num_states + 2, sigma = target_energy, which = 'LM', return_eigenvectors = False) #Get new eigenvalues
    vals = np.sort(vals)

    vals.tofile(e_list)
    e_log.write("Iter %d complete\n" % (t+1))

    e_list.flush()
    e_log.flush()

e_list.close()
e_log.close()

我把num_itr设置为100。在num_itr循环的第9次循环中（通过e_log写入了8行来表示），程序崩溃并显示错误信息

无法扩展MemType 0: jcol 7438

追踪（最近的调用最后）：

   File "/usr/lusers/clb37/QREM_Energy_Gatherer.py", line 55, in <module>

           vals = sparse.linalg.eigsh(H, k = num_states + 2, sigma = target_energy, which = 'LM', return_eigenvectors = False)

   File "/usr/lusers/clb37/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.py", line 1524, in eigsh

           symmetric=True, tol=tol)

   File "/usr/lusers/clb37/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.py", line 1030, in get_OPinv_matvec

           return SpLuInv(A.tocsc()).matvec

   File "/usr/lusers/clb37/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.py", line 898, in __init__

           self.M_lu = splu(M)

   File "/usr/lusers/clb37/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/scipy/sparse/linalg/dsolve/linsolve.py", line 242, in splu

           ilu=False, options=_options)

内存错误

确实，每次我在我的机器上运行时，程序都会在第9次循环时崩溃，当我尝试在内存更大的机器上运行这段代码时，程序可以运行更多的迭代才崩溃，所以看起来电脑确实是内存不够了。如果问题就仅仅是这样那也没关系，但我不明白的是，为什么程序在第1次迭代时不会崩溃。我没有看到num_itr循环中的8行代码有任何地方会写入内存而不在下一次迭代中被覆盖。我使用Heapy的heap()函数查看我的内存使用情况，每次都只打印出“总大小 = 11715240字节”。

我觉得这里一定有我不知道的基本问题，可能是我写代码时的某个bug，或者是关于内存处理的某个细节。有没有人能解释一下，为什么这段代码在num_itr循环的第9次时会失败，而在第1次时却不会？

Answer 1

好的，这个问题在Scipy 0.14.0版本上似乎是可以重复出现的。

看起来可以通过在循环中添加

import gc; gc.collect()

来解决这个问题，这样可以强制Python的循环垃圾回收器运行。

问题的根源在于，scipy.sparse.eigh内部某个地方存在一个循环引用的情况，类似于：

class Foo(object):
    pass
a = Foo()
b = Foo()
a.spam = b
b.spam = a
del a, b   # <- but a, b still refer to each other and are not dead

从原则上来说，这种情况是可以接受的：虽然Python的引用计数无法检测到这种循环垃圾，但系统会定期运行一次回收，来清理这些对象。不过，如果每个对象在内存中占用的空间非常大（比如大的Numpy数组），那么这种定期回收就显得太少了，可能在下一次回收之前就会耗尽内存。

因此，一个解决办法是，当你知道有大量垃圾需要清理时，强制垃圾回收器运行。更好的解决办法是修改scipy.sparse.eigh，使得一开始就不会产生这种循环垃圾。