有python库用于非标准代数类对象的稀疏矩阵操作吗？

Question

摘要：我正在寻找一种方法来使用稀疏矩阵进行计算，稀疏矩阵的非零项不是通常的整数/浮点等，而是代数的元素，即带有加法、乘法和零元素的非标准python类的实例。你知道吗

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

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

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

它不适用于稀疏矩阵。 为了加快计算速度，我现在想用稀疏矩阵代替这些密集矩阵。我曾尝试使用SciPy的2-D稀疏矩阵包scipy.sparse来实现这一点，但到目前为止我失败了。我可以用我的代数元素填充这些稀疏矩阵类的实例，但是每当我用它们进行计算时，都会收到如下错误消息

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

对我来说，这表明对scipy.sparse支持的对象类型有限制。我不认为稀疏矩阵的运算应该关心对象类型的数学原因。只要类拥有所有float操作，比如说，它就应该工作。我错过了什么？有没有一种替代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上使用Python3.5.2。你知道吗

Answer 1

这可能更多地属于评论类，但作为一个答案，我可以把它做得更长，编辑得更多。你知道吗

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

scipy.sparse还没有开发这种功能。一个coo格式的矩阵可能可以用对象输入来创建，但这是因为它做的不多。事实上，如果data、row和col输入具有正确的numpy数组设置，则它们将作为coo属性使用，而不会发生更改。你知道吗

显然，像使用csr等一样使用indptr也只是分配属性。从coo到csr的转换可能不太好，因为这涉及到重复项的求和。你知道吗

在任何情况下，csr数学代码混合使用python和c（cython），编译部分使用有限数量的数字类型-长整数、双整数和浮点。我认为它甚至不适用于短int（int8，int16）。它不实现ndarrays所做的任何对象数据类型委托。你知道吗

用你的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不需要更改值。但是回到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有自己的稀疏系列/数据帧实现

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

By default when converting to CSR or CSC format, duplicate (i,j) entries will be summed together. This facilitates efficient construction of finite element matrices and the like. (see example)