Разреженная 3d-матрица/массив в Python?

Ответ 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))

Вы можете сделать эту реализацию более надежной, проверив, что вход на самом деле является кортежем, и что он содержит только целые числа, но это просто замедлит работу, поэтому я бы не стал беспокоиться, если вы не отпустили свой код до мир позже.

EDIT - реализация 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

Несмотря на то, что вам всегда нужно вручную отсканировать это для проблемы, потому что (как указано в комментарии к коду) вам нужно будет определить, какие индексы вы суммируете, и будьте осторожны относительно длины массива или вещей, выигранных работа!

РЕДАКТИРОВАТЬ 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 для повышения производительности, прежде чем беспокоиться об этом.

Ответ 2

Посмотрите sparray - разреженные n-мерные массивы в Python (автор Jan Erik Solem). Также доступен на github.

Ответ 3

Чем лучше писать все с нуля, тем лучше, возможно, использовать scipy разреженный модуль. Это может привести к (большей) эффективности. У меня была несколько схожая проблема, но мне приходилось только эффективно обращаться к данным, а не выполнять какие-либо операции над ними. Кроме того, мои данные были разрежены только в двух из трех измерений.

Я написал класс, который решает мою проблему и может (насколько я думаю) легко расширить, чтобы удовлетворить потребности 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])])

Ответ 4

Альтернативным ответом на этот год является 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>