是否有用于非标准代数类对象的稀疏矩阵运算的 python 库？

Question

总结：我正在寻找一种方法来使用稀疏矩阵进行计算，其非零条目不是通常的整数/浮点数/等，而是代数的元素，即具有加法、乘法和零元素的非标准 python 类。

它适用于密集矩阵。我已经通过定义一个 python 类 algebra 并重载加法和乘法来实现这个代数：

class algebra(object):
   ...
   __mul__(self,other):
      ...
   __add__(self,other):
      ...

numpy 允许我定义向量和矩阵，其条目是类algebra 的实例。它还允许我执行所有常用操作，如矩阵乘法/加法/张量点/切片/等，因此它的工作原理与整数/浮点数/等上的矩阵一样。

它不适用于稀疏矩阵。 为了加快计算速度，我现在想用稀疏矩阵替换这些密集矩阵。我试图用 SciPy 的二维稀疏矩阵包scipy.sparse 来完成这项工作，但到目前为止我失败了。我可以通过我的代数元素填充这些稀疏矩阵类的实例，但是每当我对它们进行计算时，我都会收到一条错误消息，例如

TypeError: no supported conversion for types: (dtype('O'),dtype('O'))

对我来说，这表明scipy.sparse 支持的对象类型存在限制。我看不出为什么稀疏矩阵的运算应该关心对象类型的任何数学原因。比如说，只要该类具有浮点数的所有操作，它就应该可以工作。我错过了什么？有没有支持任意对象类型的scipy.sparse 的替代品？

---

下面是一个最小的工作示例。请注意，我已经用通常的整数 0 实现了代数的零元素。还请注意，我感兴趣的实际代数比实整数更复杂！

import numpy as np
from scipy.sparse import csr_matrix

class algebra(object): # the algebra of the real integers

    def __init__(self,num):
        self.num = num

    def __add__(self,other):
        if isinstance(other, self.__class__):
            return algebra(self.num+other.num)
        else:
            return self

    def __radd__(self,other):
        if isinstance(other, self.__class__):
            return algebra(self.num+other.num)
        else:
            return self

    def __mul__(self,other):
        if isinstance(other, self.__class__):
            return algebra(self.num*other.num)
        else:
            return 0

    def __rmul__(self,other):
        if isinstance(other, self.__class__):
            return algebra(self.num*other.num)
        else:
            return 0

    def __repr__(self):
        return "algebra:"+str(self.num)  

a=algebra(5)
print(a*a)
print(a*0)
print(0*a)
indptr = np.array([0, 2, 3, 6])
indices = np.array([0, 2, 2, 0, 1, 2])
data = np.array([a,a,a,a,a,a])
S = csr_matrix((data, indices, indptr), shape=(3, 3))
print(S)
print("Everything works fine up to here.")
S*S    

输出是：

algebra:25
0
0
  (0, 0)    algebra:5
  (0, 2)    algebra:5
  (1, 2)    algebra:5
  (2, 0)    algebra:5
  (2, 1)    algebra:5
  (2, 2)    algebra:5
Everything works fine up to here.
Traceback (most recent call last):
  File "test", line 46, in <module>
    S*S    
  File "/usr/lib/python3/dist-packages/scipy/sparse/base.py", line 319, in __mul__
    return self._mul_sparse_matrix(other)
  File "/usr/lib/python3/dist-packages/scipy/sparse/compressed.py", line 499, in _mul_sparse_matrix
    data = np.empty(nnz, dtype=upcast(self.dtype, other.dtype))
  File "/usr/lib/python3/dist-packages/scipy/sparse/sputils.py", line 57, in upcast
    raise TypeError('no supported conversion for types: %r' % (args,))
TypeError: no supported conversion for types: (dtype('O'), dtype('O'))

我在 Linux 上使用 Python 3.5.2。

Answer 1

这可能在注释类别中变得更多，但作为一个答案，我可以更长的答案，并更多地编辑它。

numpy阵列通过将指针/引用存储对数组的数据缓冲区中的对象来实现对象dtype。数学是通过将任务委托给对象方法来完成的。迭代基本上是Python速度，与列表理解相当（可能是有点慢）。 numpy不在这些对象上进行快速编译的数学。

scipy.sparse @尚未开发出这种功能。 A coo格式矩阵可能会使用对象输入创建 - 但这是因为它没有多少。事实上，如果data和col和col输入有numpy array设置，它们用作coo属性而无需更改。

显然是在您使用indptr @等@ et等时分配属性。 A coo至csr转换可能无法工作得不好，因为涉及重复的求和。

在任何情况下，csr math代码使用python和c（cython）的混合，并且编译的部分适用于有限数量的数字类型 - long和double整数和floats。我认为它甚至不适用于短期（@ 987654338，int16）。它不实现任何对象DType委托ndarrays do。

与您的S：

In [187]: S.A                                                                                                
...
ValueError: unsupported data types in input

In [188]: S.tocoo()                                                                                          
Out[188]: 
<3x3 sparse matrix of type '<class 'numpy.object_'>'
    with 6 stored elements in COOrdinate format>

tocoo no value更改。但回到csr需要求和重复：

In [189]: S.tocoo().tocsr()                                                                                  
 ...
TypeError: no supported conversion for types: (dtype('O'),)

In [190]: S.tolil()                                                                                          
/usr/local/lib/python3.6/dist-packages/scipy/sparse/sputils.py:115: UserWarning: object dtype is not supported by sparse matrices
  warnings.warn("object dtype is not supported by sparse matrices")
Out[190]: 
<3x3 sparse matrix of type '<class 'numpy.object_'>'
    with 6 stored elements in LInked List format>

在存储此对象数据时没有问题

数学与您的对象列表与数组 - 类似的时间：

In [192]: alist = [a]*100                                                                                    
In [193]: arr = np.array(alist)                                                                              
In [194]: timeit [i*j for i,j in zip(alist,alist)]                                                           
77.9 µs ± 272 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [195]: timeit arr*arr                                                                                     
75.1 µs ± 2.29 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

早期的问题，您可能已经看到（我刚收到升空），在稀疏矩阵中使用int16。相同的基本问题：

Why can't I assign data to part of sparse matrix in the first "try:"?

符号库具有稀疏矩阵模块：https://docs.sympy.org/latest/modules/matrices/sparse.html

Pandas有自己的稀疏系列/ dataframe实现

https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.coo_matrix.html#scipy.sparse.coo_matrix

默认情况下，当转换为 CSR 或 CSC 格式时，重复的 (i,j) 条目将被汇总在一起。这有助于有效构建有限元矩阵等。 （见例子）