разреженная трехмерная матрица / массив в Python?

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

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

РЕДАКТИРОВАТЬ 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. Согласно самому пакету, он реализует разреженные многомерные массивы поверх 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>

Взгляните на sparray - разреженные n-мерные массивы в Python (Ян Эрик Солем). Также доступно на github.

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

Мне также нужна разреженная трехмерная матрица для решения двумерных уравнений теплопроводности (2 пространственных измерения плотные, но временное измерение - диагональ плюс и минус одна недиагональ.) Я нашел эта ссылка, которая поможет мне. Уловка состоит в том, чтобы создать массив Number, который отображает двумерную разреженную матрицу на одномерный линейный вектор. Затем постройте 2D-матрицу, построив список данных и индексов. Позже матрица Number используется для упорядочивания ответа обратно в двумерный массив.

[edit] После моего первого сообщения мне пришло в голову, что с этим можно было бы лучше справиться, используя метод .reshape(-1). После исследования метод reshape лучше, чем flatten, потому что он возвращает новое представление в исходный массив, но flatten копирует массив. В коде используется исходный массив Number. Я постараюсь обновить позже. [конец редактирования]

Я тестирую это, создавая одномерный случайный вектор и решая второй вектор. Затем умножьте его на разреженную 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}")