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

Posted 2023-03-12

【中文标题】Python中的稀疏3d矩阵/数组？

【英文标题】：sparse 3d matrix/array in Python?

【发布时间】：2011-12-02 20:42:41

【问题描述】：

在 scipy 中，我们可以使用 scipy.sparse.lil_matrix() 等构造一个稀疏矩阵。但是矩阵是二维的。

我想知道 Python 中是否存在用于稀疏 3d 矩阵/数组（张量）的数据结构？

附言我在 3d 中有很多稀疏数据，需要一个张量来存储/执行乘法。如果没有现有的数据结构，有什么建议可以实现这样的张量？

【问题讨论】：

这篇文章可能对***.com/questions/4490961/…有帮助

...但不是稀疏的，不幸的是。

“二维矩阵”是什么意思？如果您的意思是表示 2D 线性变换的矩阵，那么您说的是 2x2 实数值矩阵（由浮点值近似），行列式为 1，用于刚性旋转。如果您还想表示平移，那么您将 2x2 矩阵嵌入到 3x3 矩阵中，如果您想允许剪切或扩展/收缩，那么您可以放宽行列式要求 - 但即便如此，总共有 9 个浮点值。为什么想要/需要稀疏表示？

@Peter “二维矩阵”表示二维矩阵。二维矩阵中的单位可以表示为 (x,y, r)，其中 x & y 是坐标，r 是存储在 (x, y) 处的值。我需要一个稀疏表示，因为当 x & y 非常大时，比如 x

谢谢 - 我被 P.S. 弄糊涂了。在您的问题中（在我看来，您想将一堆欧几里得元组乘以矩阵，线性代数风格）。但是，如果您谈论的是 m x n x o 矩阵，那么听起来您的“稀疏”实现需要提供某种迭代器接口才能实现（逐个元素）乘法。

【参考方案1】：

很高兴提出一个（可能是显而易见的）实现，如果您有时间和空间来构建新的依赖项，并且需要它更快，则可以使用纯 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))

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

编辑 - 矩阵乘法问题的 Cython 实现，假设其他张量是 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 实现的性能。

【讨论】：

问题在于数学运算，而不是数据容器...我从未听说过高效稀疏 N-d 张量积的算法。看看scipy.sparse.dok_matrix。这就是您在此处描述的内容，仅限于 2D。扩展它来保存 N-D 数据很容易，但是如何对数据进行操作？ （话虽如此，你的回答是完全合理的......）

啊，我误会了？所以这个问题是在询问更多关于 scipy 兼容矩阵乘法运算的实现？当然，这应该相对容易实现，因为您真正需要的只是对索引处的值进行循环查询，这是我提供的。不过，我会看看 scipy 规范。

好吧，可以说，我也误解了。无论哪种方式，我的观点是您在进行操作时没有利用稀疏结构。您在编辑中描述的内容将其视为密集数组。 （这当然有效！您的答案解决了手头的问题。）稀疏矩阵库利用了数组的稀疏性，并避免了诸如遍历数组的每个元素之类的事情，而不管“稀疏性”如何。这就是使用稀疏矩阵的要点。操作大致随着“密集”元素的数量而不是矩阵的整体尺寸而缩放。

@tehwalrus 感谢您的回复。但是恐怕与您建议的数据结构相乘可能不是很有效...

@JoeKington 你必须遍历非稀疏数组中的每个元素（在这种情况下为u），当然？除非两者都是稀疏的，否则我误解得更多。在这种情况下，您可以简单地遍历字典中的键，并从元组中提取索引。无论如何，我跟不上稀疏代数的速度，更不用说优化该主题算法背后的计算机科学了。对不起，@zhongqi！

【参考方案2】：

看看sparray - sparse n-dimensional arrays in Python（作者 Jan Erik Solem）。也可通过github 获得。

【参考方案3】：

比从头开始编写所有新东西更好的方法可能是尽可能使用 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])])

【讨论】：

我也有同样的情况，这非常有用。我认为通过一些额外的工作，这可以在scipy 的稀疏数组模块中实现。你考虑过吗？

@TomCho 谢谢！我没有考虑将它实现到 scipy 的稀疏模块中。我认为 scipy 中的实现应该支持所有标准的 numpy 矩阵功能。这将是可行的，但需要大量的工作。另外，我认为为这些矩阵上的操作添加 C 实现会更有效，更适合 scipy。

【参考方案4】：

截至 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>

【讨论】：

@JayShin 也许，但老实说，我认为你必须测试一下。

我建议写“截至 2017 年”而不是“截至今年”

有什么理由指向这个分叉而不是原来的（github.com/pydata/sparse 分叉源）之一？此处链接的链接自 2018 年以来未更新，而我链接的链接最近更新（2021 年 1 月 4 日）。

@0xc0de 原因是我不够重视！我用原始链接修复了它。谢谢

【参考方案5】：

我还需要 3D 稀疏矩阵来求解 2D 热方程（2 个空间维度是密集的，但时间维度是对角线正负一非对角线。）我找到了 this 链接来指导我。诀窍是创建一个数组Number，将二维稀疏矩阵映射到一维线性向量。然后通过构建数据和索引列表来构建 2D 矩阵。稍后，Number 矩阵用于将答案排列回二维数组。

[edit] 在我第一次发帖后我突然想到，使用.reshape(-1) 方法可以更好地处理这个问题。经过研究，reshape 方法比flatten 更好，因为它返回一个新视图到原始数组中，但flatten 复制了数组。该代码使用原始的Number 数组。我稍后会尝试更新。[结束编辑]

我通过创建一维随机向量并求解第二个向量来对其进行测试。然后将它乘以稀疏二维矩阵，得到相同的结果。

注意：我在一个循环中重复多次，使用完全相同的矩阵M，所以你可能认为解决inverse(会更有效M)。但是 M 的倒数是 not 稀疏的，所以我认为（但尚未测试）使用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 = \nNumber")

    # 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 = \nMySpNodes.toarray()")
b = rand(nX*nY)
print(f"b=\nb")
x = spLA.spsolve(MySpNodes.tocsr(), b)
print(f"x=\nx")
print(f"Multiply back together=\nx * MySpNodes")

【参考方案6】：

我需要一个 x、y、z 的 3d 查找表并想出了这个解决方案。 为什么不使用其中一个维度作为第三维度的除数？ IE。使用 x 和 'yz' 作为矩阵维度

例如。如果 x 有 80 个潜在成员，y 有 100 个潜在成员，z 有 20 个潜在成员 您使稀疏矩阵到 2000 年为 80（即 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
====================