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
#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
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;
import scipy.sparse as sp
import numpy as np
class Sparse3D():
"""
Class to store and access 3 dimensional sparse matrices efficiently
"""
def __init__(self, *sparseMatrices):
"""
Constructor
Takes a stack of sparse 2D matrices with the same dimensions
"""
self.data = sp.vstack(sparseMatrices, "dok")
self.shape = (len(sparseMatrices), *sparseMatrices[0].shape)
self._dim1_jump = np.arange(0, self.shape[1]*self.shape[0], self.shape[1])
self._dim1 = np.arange(self.shape[0])
self._dim2 = np.arange(self.shape[1])
def __getitem__(self, pos):
if not type(pos) == tuple:
if not hasattr(pos, "__iter__") and not type(pos) == slice:
return self.data[self._dim1_jump[pos] + self._dim2]
else:
return Sparse3D(*(self[self._dim1[i]] for i in self._dim1[pos]))
elif len(pos) > 3:
raise IndexError("too many indices for array")
else:
if (not hasattr(pos[0], "__iter__") and not type(pos[0]) == slice or
not hasattr(pos[1], "__iter__") and not type(pos[1]) == slice):
if len(pos) == 2:
result = self.data[self._dim1_jump[pos[0]] + self._dim2[pos[1]]]
else:
result = self.data[self._dim1_jump[pos[0]] + self._dim2[pos[1]], pos[2]].T
if hasattr(pos[2], "__iter__") or type(pos[2]) == slice:
result = result.T
return result
else:
if len(pos) == 2:
return Sparse3D(*(self[i, self._dim2[pos[1]]] for i in self._dim1[pos[0]]))
else:
if not hasattr(pos[2], "__iter__") and not type(pos[2]) == slice:
return sp.vstack([self[self._dim1[pos[0]], i, pos[2]]
for i in self._dim2[pos[1]]]).T
else:
return Sparse3D(*(self[i, self._dim2[pos[1]], pos[2]]
for i in self._dim1[pos[0]]))
def toarray(self):
return np.array([self[i].toarray() for i in range(self.shape[0])])
4 回答
很高兴建议一个（可能是显而易见的）实现，如果你有新的依赖项的时间和空间，可以在纯Python或C / Cython中进行，并且需要它更快 .
N维中的稀疏矩阵可以假设大多数元素都是空的，因此我们使用键入元组的字典：
你会像这样使用它：
您可以通过验证输入实际上是一个元组，并且它只包含整数来使这个实现更加健壮，但这只会减慢速度，所以除非您以后将代码发布到全世界，否则我不会担心 .
EDIT - 矩阵乘法问题的Cython实现，假设其他张量是N维NumPy数组（
numpy.ndarray）可能如下所示：虽然你总是需要手动解决这个问题，因为（如代码注释中所述）你需要定义你正在总结的索引，并且要小心数组长度或事情不起作用！
EDIT 2 - 如果另一个矩阵也是稀疏的，那么您不需要进行三向循环：
我对C实现的建议是使用一个简单的结构来保存索引和值：
然后，您需要一些函数来分配和维护这些结构的动态数组，并根据需要快速搜索它们;但是你应该在担心这些东西之前测试Python实现的性能 .
看看sparray - sparse n-dimensional arrays in Python（由Jan Erik Solem撰写） . 也可在github上找到 .
截至今年的另一个答案是sparse包 . 根据包本身，它通过推广
scipy.sparse.coo_matrix布局在NumPy和scipy.sparse之上实现稀疏多维数组 .以下是从文档中获取的示例：
比从头开始写新东西更好，可能是尽可能使用scipy的稀疏模块 . 这可能会导致（更多）更好的性能 . 我有一个类似的问题，但我只需要有效地访问数据，而不是对它们执行任何操作 . 此外，我的数据在三个维度中只有两个是稀疏的 .
我写了一个解决我问题的课程，可以（据我认为）轻松扩展以满足OP的需求 . 不过，它可能仍然有一些改进的潜力 .