matriz - sparse matrix python

matriz/matriz 3D dispersa en Python?

Question

En scipy, podemos construir una matriz dispersa usando scipy.sparse.lil_matrix () etc. Pero la matriz está en 2d.

Me pregunto si existe una estructura de datos existente para matrix / array escasa (tensor) en Python.

ps Tengo muchos datos escasos en 3D y necesito un tensor para almacenar / realizar la multiplicación. ¿Alguna sugerencia para implementar dicho tensor si no hay una estructura de datos existente?

Answer 1

Eche un vistazo a las matrices n-dimensionales dispersas en sparray en Python (por Jan Erik Solem). También disponible en github .

Answer 2

Feliz de sugerir una implementación (posiblemente obvia) de esto, que podría hacerse en Python puro o C / Cython si tiene tiempo y espacio para nuevas dependencias, y necesita que sea más rápido.

Una matriz dispersa en N dimensiones puede suponer que la mayoría de los elementos están vacíos, por lo que usamos un diccionario con clave en tuplas:

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

y lo usarías así:

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

Podrías hacer que esta implementación sea más robusta al verificar que la entrada es de hecho una tupla, y que contiene solo enteros, pero que ralentizará las cosas para que no me preocupe a menos que liberes tu código al mundo más tarde.

EDITAR - una implementación de Cython del problema de multiplicación de la matriz, suponiendo que otro tensor es una matriz numpy.ndarray NumPy ( numpy.ndarray ) podría verse así:

#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

Aunque siempre tendrás que pasar este problema manualmente, porque (como se menciona en el comentario del código) necesitarás definir qué índices estás sumando, y tener cuidado con las longitudes de la matriz o las cosas no funcionarán. !

EDIT 2 - si la otra matriz también es escasa, entonces no es necesario hacer el bucle de tres vías:

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

Mi sugerencia para una implementación de C sería usar una estructura simple para contener los índices y los valores:

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

luego necesitará algunas funciones para asignar y mantener una matriz dinámica de tales estructuras, y buscarlas tan rápido como lo necesite; pero debes probar la implementación de Python para el rendimiento antes de preocuparte por eso.

Answer 3

Más lindo que escribir todo lo nuevo desde cero puede ser utilizar el módulo disperso de scipy lo más lejos posible. Esto puede conducir a (mucho) mejor rendimiento. Tuve un problema similar, pero solo tuve que acceder a los datos de manera eficiente, no realizar ninguna operación con ellos. Además, mis datos solo eran escasos en dos de las tres dimensiones.

He escrito una clase que resuelve mi problema y podría (hasta donde yo creo) extenderse fácilmente para satisfacer las necesidades del OP. Sin embargo, aún puede tener algún potencial de mejora.

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 4

Una respuesta alternativa a partir de este año es el paquete sparse . De acuerdo con el paquete en sí, implementa arrays multidimensionales dispersos sobre NumPy y scipy.sparse al generalizar el diseño de scipy.sparse.coo_matrix .

Aquí hay un ejemplo tomado de los documentos:

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>