我很乐意提出一个（可能是显而易见的）实现，如果你有时间和空间来处理新的依赖项，并且需要它更快，可以用纯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))

型
您可以通过验证输入实际上是一个元组，并且它只包含整数来使这个实现更健壮，但这只会减慢速度，所以我不会担心，除非您稍后将代码发布给世界。

EDIT-矩阵乘法问题的Cython实现，假设其他Tensor是N维NumPy数组（numpy.ndarray）可能看起来像这样：

#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

型
尽管你总是需要手动处理手头的问题，因为（如代码注解中所提到的）你需要定义你在哪些索引上求和，并且要小心数组长度，否则事情就不起作用了！

EDIT 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实现的性能。

截至2017年的另一个答案是sparse包。根据包本身，它通过泛化scipy.sparse.coo_matrix布局在NumPy和scipy.sparse之上实现稀疏多维数组。
下面是一个从文档中提取的示例：

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>

字符串

看看sparray - sparse n-dimensional arrays in Python（作者：Jan Erik Solem）。也可以在github上找到。

比从头开始编写所有新内容更好的方法是尽可能使用scipy的sparse模块。这可能会带来（更）更好的性能。我遇到了一个有点类似的问题，但我只需要有效地访问数据，而不是对它们执行任何操作。此外，我的数据只在三维中的二维中稀疏。
我已经写了一个类来解决我的问题，并且（据我所知）可以很容易地扩展以满足OP的需要。

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])])

字符串

我还需要3D稀疏矩阵求解2D热方程（2个空间维度是密集的，但时间维度是对角线加和减一个非对角线。我找到了this链接来指导我。诀窍是创建一个数组Number，将2D稀疏矩阵Map到1D线性向量。然后通过构建数据和索引列表来构建2D矩阵。稍后，Number矩阵用于将答案排列回2D数组。
[edit]在我第一次发表文章后，我突然想到，使用.reshape(-1)方法可以更好地处理这个问题。经过研究，reshape方法比flatten方法更好，因为它返回一个新的视图到原始数组中，但是flatten复制了数组。代码使用原始Number数组。我稍后会尝试更新。[end edit]
我通过创建一个1D随机向量并求解第二个向量来测试它，然后将其乘以稀疏的2D矩阵，我得到了相同的结果。

备注：我在一个循环中用完全相同的矩阵M重复了很多次，所以你可能会认为求解inverse(M)会更有效。但是M的逆矩阵 * 不是 * 稀疏的，所以我认为（但没有测试过）使用spsolve是一个更好的解决方案。“最好”可能取决于你使用的矩阵有多大。

#!/usr/bin/env python3
# testSparse.py
# profhuster

import numpy as np
import scipy.sparse as sM
import scipy.sparse.linalg as spLA
from array import array
from numpy.random import rand, seed
seed(101520)

nX = 4
nY = 3
r = 0.1

def loadSpNodes(nX, nY, r):
    # Matrix to map 2D array of nodes to 1D array
    Number = np.zeros((nY, nX), dtype=int)

    # Map each element of the 2D array to a 1D array
    iM = 0
    for i in range(nX):
        for j in range(nY):
            Number[j, i] = iM
            iM += 1
    print(f"Number = \n{Number}")

    # Now create a sparse matrix of the "stencil"
    diagVal = 1 + 4 * r
    offVal = -r
    d_list = array('f')
    i_list = array('i')
    j_list = array('i')
    # Loop over the 2D nodes matrix
    for i in range(nX):
        for j in range(nY):
            # Recall the 1D number
            iSparse = Number[j, i]
            # populate the diagonal
            d_list.append(diagVal)
            i_list.append(iSparse)
            j_list.append(iSparse)
            # Now, for each rectangular neighbor, add the 
            # off-diagonal entries
            # Use a try-except, so boundry nodes work
            for (jj,ii) in ((j+1,i),(j-1,i),(j,i+1),(j,i-1)):
                try:
                    iNeigh = Number[jj, ii]
                    if jj >= 0 and ii >=0:
                        d_list.append(offVal)
                        i_list.append(iSparse)
                        j_list.append(iNeigh)
                except IndexError:
                    pass
    spNodes = sM.coo_matrix((d_list, (i_list, j_list)), shape=(nX*nY,nX*nY))
    return spNodes

MySpNodes = loadSpNodes(nX, nY, r)
print(f"Sparse Nodes = \n{MySpNodes.toarray()}")
b = rand(nX*nY)
print(f"b=\n{b}")
x = spLA.spsolve(MySpNodes.tocsr(), b)
print(f"x=\n{x}")
print(f"Multiply back together=\n{x * MySpNodes}")

字符串

我需要一个三维查找表的x，y，z和想出了这个解决方案。
为什么不使用其中一个维度作为第三个维度的除数呢？即使用x和'yz'作为矩阵维度
例如，如果x有80个潜在成员，y有100个潜在成员，z有20个潜在成员，那么稀疏矩阵为80 × 2000（即xy=100x20）
x维度与往常一样
yz维度：前100个元素将表示z=0，y=0到99
...第二个100将表示z=2，y=0到99等
所以位于（x，y，z）给定元素将在（x，z*100 + y）处的稀疏矩阵中
如果你需要使用负数，设计一个任意的偏移量到你的矩阵转换中。如果需要，这个解决方案可以扩展到n维

from scipy import sparse
m = sparse.lil_matrix((100,2000), dtype=float)

def add_element((x,y,z), element):
    element=float(element)
    m[x,y+z*100]=element

def get_element(x,y,z):
    return m[x,y+z*100]

add_element([3,2,4],2.2)
add_element([20,15,7], 1.2)
print get_element(0,0,0)
print get_element(3,2,4)
print get_element(20,15,7)
print "  This is m sparse:";print m

====================
OUTPUT:
0.0
2.2
1.2
  This is m sparse:
  (3, 402L) 2.2
  (20, 715L)    1.2
====================

字符串

np.nditer也可以（错误）用于创建数组的DOC风格的稀疏表示：
在我的用例中，我需要一个所有单元格!= 0的稀疏数组。

import numpy as np

def dok_ndimensional_sparse(nd_histogram: np.ndarray) -> dict[tuple, np.ndarray]:
    """DOK (Dictionary of Keys) style sparse array of n-dimsional array."""
    with np.nditer(nd_histogram, flags=["multi_index"]) as it:
        return {it.multi_index: cell_value for cell_value in it if cell_value != 0}

x = np.array([[[0., 0.],
               [0., 1.]]])

dok_ndimensional_sparse(x)
# {(0, 1, 1): array(1.)}

字符串