Pickling scipy 的 SuperLU 类用于不完全 LU 分解答案

【问题标题】：Pickling scipy's SuperLU class for incomplete LU factorizationPickling scipy 的 SuperLU 类用于不完全 LU 分解
【发布时间】：2015-06-19 16:15:34
【问题描述】：

使用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 分解的结果缓存到文件中吗？

【问题讨论】：

不幸的是，我没有看到任何直接的方法来腌制SuperLU 实例。如果您可以直接实例化SuperLU，那么您可以使用copy_reg.pickle 来注册自定义酸洗函数，但正如您所发现的，这似乎是不可能的。正如@MikeMcKerns 提到的，您可以保存L 和U 稀疏三角矩阵，但为了使用它们，您仍然需要自己重新实现SuperLU.solve()。
这听起来确实是一个不错的功能，所以我建议您在 the scipy issue tracker 上发布新功能请求
好点@ali_m，我刚刚提交了一个请求。

标签： python numpy scipy linear-algebra pickle

【解决方案1】：

作为实现.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]

【讨论】：

【解决方案2】：

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

【讨论】：

我有同样的问题/需要将稀疏 lu 因子的结果保存到磁盘并恢复。我看到你 3 年前的功能请求仍然开放。你有没有找到其他直接访问求解函数的方法？

【解决方案3】：

如果你能找到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)
'\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 的内置序列化（泡菜更小），但我不知道。

【讨论】：

不幸的是SuperLU 类确实是implemented purely in C，所以酸洗是一个非首发。正如 OP 所提到的，SuperLU 实例也不能直接实例化，因此除非您自己重新实现 .solve() 方法，否则您对下三角稀疏矩阵和上三角稀疏矩阵无能为力。