matrice/réseau 3d clairsemé en Python ?

Question

En scipy, nous pouvons construire une matrice éparse en utilisant scipy.sparse.lil_matrix() etc. Mais la matrice est en 2d.

Je me demande s'il existe une structure de données pour les matrices / tableaux 3d éparses (tenseurs) en Python ?

p.s. J'ai beaucoup de données éparses en 3D et j'ai besoin d'un tenseur pour stocker / effectuer des multiplications. Avez-vous des suggestions pour implémenter un tel tenseur s'il n'y a pas de structure de données existante ?

Answer 1

Je serais heureux de suggérer une implémentation (peut-être évidente) de ceci, qui pourrait être faite en Python pur ou en C/Cython si vous avez le temps et l'espace pour de nouvelles dépendances, et si vous avez besoin que ce soit plus rapide.

Une matrice éparse en N dimensions peut supposer que la plupart des éléments sont vides, nous utilisons donc un dictionnaire avec des clés sur les tuples :

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

et vous l'utiliserez comme ça :

sparse = NDSparseMatrix()
sparse.addValue((1,2,3), 15.7)
should_be_zero = sparse.readValue((1,5,13))

Vous pourriez rendre cette implémentation plus robuste en vérifiant que l'entrée est en fait un tuple, et qu'elle ne contient que des entiers, mais cela ne fera que ralentir les choses, donc je ne m'en préoccuperais pas à moins que vous ne diffusiez votre code au monde entier plus tard.

EDIT - une implémentation Cython du problème de multiplication de matrice, en supposant que l'autre tenseur est un tableau NumPy à N dimensions ( numpy.ndarray ) pourrait ressembler à ceci :

#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

Bien que vous aurez toujours besoin de le faire à la main pour le problème en question, parce que (comme mentionné dans le commentaire du code) vous aurez besoin de définir quels indices vous additionnez, et faites attention aux longueurs des tableaux ou les choses ne fonctionneront pas !

EDIT 2 - si l'autre matrice est également éparse, alors vous n'avez pas besoin de faire la boucle à trois voies :

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

Ma suggestion pour une implémentation en C serait d'utiliser un simple struct pour contenir les indices et les valeurs :

typedef struct {
  int index[3];
  float value;
} entry_t;

vous aurez alors besoin de fonctions pour allouer et maintenir un tableau dynamique de ces structs, et les rechercher aussi vite que nécessaire ; mais vous devriez tester l'implémentation Python en place pour les performances avant de vous préoccuper de ce genre de choses.

Answer 2

Une réponse alternative à partir de 2017 est le sparse paquet. Selon le paquet lui-même, il met en œuvre des tableaux multidimensionnels épars au-dessus de NumPy et de scipy.sparse en généralisant le scipy.sparse.coo_matrix disposition.

Voici un exemple tiré de la documentation :

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>

Answer 3

Jetez un coup d'œil à sparray - tableaux n-dimensionnels épars en Python (par Jan Erik Solem). Également disponible sur github .

Answer 4

Il est préférable d'utiliser le module sparse de Scipy dans la mesure du possible plutôt que d'écrire tout à partir de zéro. Cela peut conduire à de (bien) meilleures performances. J'ai eu un problème un peu similaire, mais je devais seulement accéder efficacement aux données, sans effectuer d'opérations sur elles. De plus, mes données n'étaient éparses que dans deux des trois dimensions.

J'ai écrit une classe qui résout mon problème et qui pourrait (pour autant que je le pense) être facilement étendue pour satisfaire les besoins de l'OP. Cependant, elle peut encore être améliorée.

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

Answer 5

J'ai également besoin d'une matrice 3D peu dense pour résoudre les équations de chaleur en 2D (les deux dimensions spatiales sont denses, mais la dimension temporelle est diagonale plus et moins une offdiagonale). J'ai trouvé ce lien pour me guider. L'astuce consiste à créer un tableau Number qui transforme la matrice 2D éparse en un vecteur linéaire 1D. Ensuite, on construit la matrice 2D en établissant une liste de données et d'indices. Plus tard, la Number est utilisée pour réorganiser la réponse en un tableau 2D.

[ modifier Il m'est venu à l'esprit, après mon message initial, que cela pourrait être mieux géré en utilisant la fonction .reshape(-1) méthode. Après la recherche, le reshape est plus efficace que la méthode flatten car elle renvoie une nouvelle vue dans le tableau original, mais flatten copie le tableau. Le code utilise l'original Number réseau. Je vais essayer de mettre à jour plus tard . montage final ]

Je le teste en créant un vecteur aléatoire 1D et en le résolvant pour un second vecteur. Je le multiplie ensuite par la matrice 2D éparse et j'obtiens le même résultat.

Note : Je répète ceci plusieurs fois en boucle avec exactement la même matrice. M donc vous pourriez penser qu'il serait plus efficace de résoudre pour inverse( M ) . Mais l'inverse de M est no clairsemé, donc je pense (mais je n'ai pas testé) qu'utiliser spsolve est une meilleure solution. La "meilleure" dépend probablement de la taille de la matrice que vous utilisez.

#!/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}")