Python中的稀疏3D矩阵/数组？

Question

在scipy这个库里，我们可以用scipy.sparse.lil_matrix()等方法来构建稀疏矩阵，不过这些矩阵都是二维的。

我想知道在Python里有没有现成的可以用来表示稀疏三维矩阵（也叫张量）的数据结构？

另外，我有很多三维的稀疏数据，需要一个张量来存储和进行乘法运算。如果没有现成的数据结构，有什么建议可以实现这样的张量吗？

Answer 1

看看这个链接：sparray - 在Python中使用稀疏n维数组（作者是Jan Erik Solem）。这个项目也可以在github上找到。

Answer 2

2017年有一个替代方案，就是 sparse 这个包。根据这个包的介绍，它是在NumPy和 scipy.sparse 的基础上实现了稀疏的多维数组，主要是对 scipy.sparse.coo_matrix 的布局进行了扩展。

下面是从文档中摘录的一个例子：

import numpy as np
n = 1000
ndims = 4
nnz = 1000000
coords = np.random.randint(0, n - 1, size=(ndims, nnz))
data = np.random.random(nnz)

import sparse
x = sparse.COO(coords, data, shape=((n,) * ndims))
x
# <COO: shape=(1000, 1000, 1000, 1000), dtype=float64, nnz=1000000>

x.nbytes
# 16000000

y = sparse.tensordot(x, x, axes=((3, 0), (1, 2)))

y
# <COO: shape=(1000, 1000, 1000, 1000), dtype=float64, nnz=1001588>

Answer 3

很高兴能给出一个可能比较明显的实现方法，你可以用纯Python或者C/Cython来做，如果你有时间和空间来添加新的依赖库，并且需要更快的速度。

在N维的稀疏矩阵中，大部分元素都是空的，所以我们可以用一个字典，字典的键是元组：

class NDSparseMatrix:
  def __init__(self):
    self.elements = {}

  def addValue(self, tuple, value):
    self.elements[tuple] = value

  def readValue(self, tuple):
    try:
      value = self.elements[tuple]
    except KeyError:
      # could also be 0.0 if using floats...
      value = 0
    return value

然后你可以这样使用它：

sparse = NDSparseMatrix()
sparse.addValue((1,2,3), 15.7)
should_be_zero = sparse.readValue((1,5,13))

你可以通过验证输入确实是一个元组，并且只包含整数，来让这个实现更健壮，但这样会让速度变慢，所以如果你不是要把代码发布给大家，就不用太担心这个。

编辑 - 如果假设另一个张量是一个N维的NumPy数组（numpy.ndarray），那么Cython实现的矩阵乘法问题可能看起来像这样：

#cython: boundscheck=False
#cython: wraparound=False

cimport numpy as np

def sparse_mult(object sparse, np.ndarray[double, ndim=3] u):
  cdef unsigned int i, j, k

  out = np.ndarray(shape=(u.shape[0],u.shape[1],u.shape[2]), dtype=double)

  for i in xrange(1,u.shape[0]-1):
    for j in xrange(1, u.shape[1]-1):
      for k in xrange(1, u.shape[2]-1):
        # note, here you must define your own rank-3 multiplication rule, which
        # is, in general, nontrivial, especially if LxMxN tensor...

        # loop over a dummy variable (or two) and perform some summation:
        out[i,j,k] = u[i,j,k] * sparse((i,j,k))

  return out

虽然你总是需要根据具体问题来手动调整这个实现，因为（正如代码注释中提到的）你需要定义你要对哪些索引进行求和，并且要小心数组的长度，否则就会出问题！

编辑 2 - 如果另一个矩阵也是稀疏的，那么你就不需要进行三重循环：

def sparse_mult(sparse, other_sparse):

  out = NDSparseMatrix()

  for key, value in sparse.elements.items():
    i, j, k = key
    # note, here you must define your own rank-3 multiplication rule, which
    # is, in general, nontrivial, especially if LxMxN tensor...

    # loop over a dummy variable (or two) and perform some summation 
    # (example indices shown):
    out.addValue(key) = out.readValue(key) + 
      other_sparse.readValue((i,j,k+1)) * sparse((i-3,j,k))

  return out

我建议在C语言中实现时，使用一个简单的结构体来存储索引和值：

typedef struct {
  int index[3];
  float value;
} entry_t;

然后你需要一些函数来分配和维护这样结构体的动态数组，并快速搜索它们；但在担心这些之前，你应该先测试一下Python实现的性能。