Pickling scipy 的 SuperLU class 用于不完全 LU 分解

Question

使用 scipy.sparse.linalg.spilu, I compute the incomplete LU-factorization of a very large sparse matrix. As this process is time-consuming I would like to save the computed LU factorization. The function returns a scipy.sparse.linalg.SuperLU 对象。

我的第一次尝试是使用 pickle 模块来保存整个对象。但是，我得到一个：

cPickle.PicklingError: Can't pickle <type 'SuperLU'>:
   attribute lookup __builtin__.SuperLU failed

错误信息。

我的第二个想法是保存 SuperLU 对象 ('L', 'U', 'nnz', 'perm_c', 'perm_r', 'shape') 的相关 class 成员，然后重新组合它。然而，SuperLU 对象似乎是不可实例化的：

>>> SuperLU()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: cannot create 'SuperLU' instances

有人知道如何将不完整的 LU 分解结果缓存到文件中吗？

Answer 1

如果你能找到 SuperLU class 的实际定义位置，如果可以 pickle... 除非 SuperLU class 在 中定义C——这是完全可能的。如果 class 在 python 中定义，那么它可以用 dill 模块进行 pickle。如果它在 C 中，那么直接 picking 那个对象就不走运了。

这是问题所在：

>>> import dill
>>> import scipy
>>> import scipy.sparse
>>> import scipy.sparse.linalg                       
>>> import numpy
>>> a = numpy.array([[1,2,2],[4,2,0],[2,0,0]])
>>> lu = scipy.sparse.linalg.splu(a)
>>> dill.detect.errors(lu)
PicklingError("Can't pickle <type 'SuperLU'>: it's not found as __builtin__.SuperLU",)
>>> lu
<SuperLU object at 0x106eb0360>
>>> lu.__class__ 
<type 'SuperLU'>
>>> lu.__class__.__module__
'__builtin__'

那你是做什么的？

我不确定这是否是完整答案，但您可以只转储构成 SparseLU 实例的稀疏矩阵。我可能遗漏了一些状态，但这应该传达了这个想法……

>>> dill.dumps(lu.L)
'\x80\x02cscipy.sparse.csc\ncsc_matrix\nq\x00)\x81q\x01}q\x02(U\x06formatq\x03U\x03cscq\x04U\x06_shapeq\x05K\x03K\x03\x86q\x06U\x06indptrq\x07cnumpy.core.multiarray\n_reconstruct\nq\x08cnumpy\nndarray\nq\tK\x00\x85q\nU\x01bq\x0b\x87q\x0cRq\r(K\x01K\x04\x85q\x0ecnumpy\ndtype\nq\x0fU\x02i4q\x10K\x00K\x01\x87q\x11Rq\x12(K\x03U\x01<q\x13NNNJ\xff\xff\xff\xffJ\xff\xff\xff\xffK\x00tq\x14b\x89U\x10\x00\x00\x00\x00\x01\x00\x00\x00\x03\x00\x00\x00\x04\x00\x00\x00q\x15tq\x16bU\x07indicesq\x17h\x08h\tK\x00\x85q\x18h\x0b\x87q\x19Rq\x1a(K\x01K\x04\x85q\x1bh\x12\x89U\x10\x00\x00\x00\x00\x01\x00\x00\x00\x02\x00\x00\x00\x02\x00\x00\x00q\x1ctq\x1dbU\x08maxprintq\x1eK2U\x04dataq\x1fh\x08h\tK\x00\x85q h\x0b\x87q!Rq"(K\x01K\x04\x85q#h\x0fU\x02f8q$K\x00K\x01\x87q%Rq&(K\x03U\x01<q\'NNNJ\xff\xff\xff\xffJ\xff\xff\xff\xffK\x00tq(b\x89U \x00\x00\x00\x00\x00\x00\xf0?\x00\x00\x00\x00\x00\x00\xf0?\x00\x00\x00\x00\x00\x00\xe0?\x00\x00\x00\x00\x00\x00\xf0?q)tq*bub.'
>>> dill.dumps(lu.U)
'\x80\x02cscipy.sparse.csc\ncsc_matrix\nq\x00)\x81q\x01}q\x02(U\x06formatq\x03U\x03cscq\x04U\x06_shapeq\x05K\x03K\x03\x86q\x06U\x06indptrq\x07cnumpy.core.multiarray\n_reconstruct\nq\x08cnumpy\nndarray\nq\tK\x00\x85q\nU\x01bq\x0b\x87q\x0cRq\r(K\x01K\x04\x85q\x0ecnumpy\ndtype\nq\x0fU\x02i4q\x10K\x00K\x01\x87q\x11Rq\x12(K\x03U\x01<q\x13NNNJ\xff\xff\xff\xffJ\xff\xff\xff\xffK\x00tq\x14b\x89U\x10\x00\x00\x00\x00\x01\x00\x00\x00\x03\x00\x00\x00\x06\x00\x00\x00q\x15tq\x16bU\x07indicesq\x17h\x08h\tK\x00\x85q\x18h\x0b\x87q\x19Rq\x1a(K\x01K\x06\x85q\x1bh\x12\x89U\x18\x00\x00\x00\x00\x00\x00\x00\x00\x01\x00\x00\x00\x00\x00\x00\x00\x01\x00\x00\x00\x02\x00\x00\x00q\x1ctq\x1dbU\x08maxprintq\x1eK2U\x04dataq\x1fh\x08h\tK\x00\x85q h\x0b\x87q!Rq"(K\x01K\x06\x85q#h\x0fU\x02f8q$K\x00K\x01\x87q%Rq&(K\x03U\x01<q\'NNNJ\xff\xff\xff\xffJ\xff\xff\xff\xffK\x00tq(b\x89U0\x00\x00\x00\x00\x00\x00\x00@\x00\x00\x00\x00\x00\x00\xf0?\x00\x00\x00\x00\x00\x00\x10@\x00\x00\x00\x00\x00\x00\x00@\x00\x00\x00\x00\x00\x00\x00@\x00\x00\x00\x00\x00\x00\xf0\xbfq)tq*bub.'

显然，您想使用 dump 而不是 dumps 来腌制文件。您甚至可以使用 numpy 的内置序列化（pickle 更小），但我不知道。

Answer 2

最后我听从了 ali_m 的建议并编写了自己的 solve() 函数。重建置换矩阵 Pr 和 Pc 后，我可以将它们与 L 和 R 一起转储并拥有我需要的一切。我还填写了 scipy 的功能请求，希望在未来的版本中会有更直接的选项。

Answer 3

作为实施.solve方法的起点... 这不太正确，因为在第一次 spsolve_triangular 调用期间有关于使用 CSC 而不是 CSR 矩阵的效率警告，但如果 L 矩阵将被多次使用，也许它是有意义的在将其存储到磁盘之前继续进行转换。

def spsolve_lu(L, U, b, perm_c=None, perm_r=None):
""" an attempt to use SuperLU data to efficiently solve
    Ax = Pr.T L U Pc.T x = b
     - note that L from SuperLU is in CSC format solving for c
       results in an efficiency warning
    Pr . A . Pc = L . U
    Lc = b      - forward solve for c
     c = Ux     - then back solve for x
"""
if perm_r is not None:
    b_old = b.copy()
    for old_ndx, new_ndx in enumerate(perm_r):
        b[new_ndx] = b_old[old_ndx]
try:    # unit_diagonal is a new kw
    c = spsolve_triangular(L, b, lower=True, unit_diagonal=True)
except TypeError:
    c = spsolve_triangular(L, b, lower=True)
px = spsolve_triangular(U, c, lower=False)
if perm_c is None:
    return px
return px[perm_c]

Pickling scipy 的 SuperLU class 用于不完全 LU 分解

Pickling scipy's SuperLU class for incomplete LU factorization

python

numpy

linear-algebra

pickle

scipy