Pickling scipy 的 SuperLU 类用于不完整的 LU 分解

Question

使用scipy.sparse.linalg.spilu，我计算了一个非常大的稀疏矩阵的不完全LU分解。由于此过程非常耗时，因此我想保存计算的LU分解。该函数返回一个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')然后重新组装它。但是，SuperLU 对象似乎是不可实例化的：

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

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

Answer 1

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

Answer 2

如果你能找到实际定义SuperLU类的位置，如果可以腌制......除非SuperLU类是用 C 语言定义的——这是完全可能的。 如果类是在 python 中定义的，那么可以使用 dill 模块对其进行腌制。 如果它在 C 中，则直接腌制该对象就不走运了。

问题是这样的：

>>> 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)
'x80x02cscipy.sparse.cscncsc_matrixnqx00)x81qx01}qx02(Ux06formatqx03Ux03cscqx04Ux06_shapeqx05Kx03Kx03x86qx06Ux06indptrqx07cnumpy.core.multiarrayn_reconstructnqx08cnumpynndarraynqtKx00x85qnUx01bqx0bx87qx0cRqr(Kx01Kx04x85qx0ecnumpyndtypenqx0fUx02i4qx10Kx00Kx01x87qx11Rqx12(Kx03Ux01<qx13NNNJxffxffxffxffJxffxffxffxffKx00tqx14bx89Ux10x00x00x00x00x01x00x00x00x03x00x00x00x04x00x00x00qx15tqx16bUx07indicesqx17hx08htKx00x85qx18hx0bx87qx19Rqx1a(Kx01Kx04x85qx1bhx12x89Ux10x00x00x00x00x01x00x00x00x02x00x00x00x02x00x00x00qx1ctqx1dbUx08maxprintqx1eK2Ux04dataqx1fhx08htKx00x85q hx0bx87q!Rq"(Kx01Kx04x85q#hx0fUx02f8q$Kx00Kx01x87q%Rq&(Kx03Ux01<q'NNNJxffxffxffxffJxffxffxffxffKx00tq(bx89U x00x00x00x00x00x00xf0?x00x00x00x00x00x00xf0?x00x00x00x00x00x00xe0?x00x00x00x00x00x00xf0?q)tq*bub.'
>>> dill.dumps(lu.U)
'x80x02cscipy.sparse.cscncsc_matrixnqx00)x81qx01}qx02(Ux06formatqx03Ux03cscqx04Ux06_shapeqx05Kx03Kx03x86qx06Ux06indptrqx07cnumpy.core.multiarrayn_reconstructnqx08cnumpynndarraynqtKx00x85qnUx01bqx0bx87qx0cRqr(Kx01Kx04x85qx0ecnumpyndtypenqx0fUx02i4qx10Kx00Kx01x87qx11Rqx12(Kx03Ux01<qx13NNNJxffxffxffxffJxffxffxffxffKx00tqx14bx89Ux10x00x00x00x00x01x00x00x00x03x00x00x00x06x00x00x00qx15tqx16bUx07indicesqx17hx08htKx00x85qx18hx0bx87qx19Rqx1a(Kx01Kx06x85qx1bhx12x89Ux18x00x00x00x00x00x00x00x00x01x00x00x00x00x00x00x00x01x00x00x00x02x00x00x00qx1ctqx1dbUx08maxprintqx1eK2Ux04dataqx1fhx08htKx00x85q hx0bx87q!Rq"(Kx01Kx06x85q#hx0fUx02f8q$Kx00Kx01x87q%Rq&(Kx03Ux01<q'NNNJxffxffxffxffJxffxffxffxffKx00tq(bx89U0x00x00x00x00x00x00x00@x00x00x00x00x00x00xf0?x00x00x00x00x00x00x10@x00x00x00x00x00x00x00@x00x00x00x00x00x00x00@x00x00x00x00x00x00xf0xbfq)tq*bub.'

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

Answer 3

作为实现.solve方法的起点...这不太正确，因为在第一次调用期间有一个关于使用 CSC 而不是 CSR 矩阵的效率警告spsolve_triangular但如果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]