Numpy matriz simétrica 'inteligente'

Question

¿Hay una matriz simétrica inteligente y eficiente en el espacio en numpy que automáticamente (y transparente) ocupa la posición en [j][i] cuando se escribe [i][j]?

import numpy 
a = numpy.symmetric((3, 3)) 
a[0][1] = 1 
a[1][0] == a[0][1] 
# True 
print(a) 
# [[0 1 0], [1 0 0], [0 0 0]] 

assert numpy.all(a == a.T) # for any symmetric matrix 

Un Hermitian automático también sería bueno, aunque no necesitaré eso al momento de escribir.

Answer 1

Si usted puede permitirse el lujo de simetrizar la matriz justo antes de hacer cálculos, lo siguiente debe ser razonablemente rápido:

def symmetrize(a): 
    return a + a.T - numpy.diag(a.diagonal()) 

Esto funciona bajo supuestos razonables (como no hacer tanto a[0, 1] = 42 y lo contradictorio a[1, 0] = 123 antes de ejecutar symmetrize)

Si realmente necesita una simetrización transparente, es posible considerar la subclasificación numpy.ndarray y simplemente redefiniendo __setitem__:

class SymNDArray(numpy.ndarray): 
    def __setitem__(self, (i, j), value): 
     super(SymNDArray, self).__setitem__((i, j), value)      
     super(SymNDArray, self).__setitem__((j, i), value)      

def symarray(input_array): 
    """ 
    Returns a symmetrized version of the array-like input_array. 
    Further assignments to the array are automatically symmetrized. 
    """ 
    return symmetrize(numpy.asarray(input_array)).view(SymNDArray) 

# Example: 
a = symarray(numpy.zeros((3, 3))) 
a[0, 1] = 42 
print a # a[1, 0] == 42 too! 

(o el equivalente con matrices en lugar de matrices, dependiendo de sus necesidades). Este enfoque incluso maneja asignaciones más complicadas, como a[:, 1] = -1, que establece correctamente a[1, :] elementos.

Tenga en cuenta que Python 3 eliminado la posibilidad de escribir def …(…, (i, j),…), por lo que el código tiene que ser adaptado ligeramente antes de ejecutar con Python 3: def __setitem__(self, indexes, value): (i, j) = indexes ...

Answer 2

La cuestión más general de tratamiento óptimo de matrices simétricas en numpy me molestó demasiado .

Después de estudiarlo, creo que la respuesta es que ese numpy está algo restringido por la disposición de memoria soportada por las rutinas BLAS subyacentes para las matrices simétricas.

Mientras que algunas rutinas BLAS explotan la simetría para acelerar los cálculos en matrices simétricas, todavía usan la misma estructura de memoria que una matriz completa, es decir, n^2 espacio en lugar de n(n+1)/2. Solo les dicen que la matriz es simétrica y que solo usan los valores en el triángulo superior o inferior.

Algunos de los scipy.linalg rutinas ACEPTA banderas (como sym_pos=True en linalg.solve), que sea traspasada a las rutinas BLAS, aunque más apoyo para esto en numpy sería bueno, en envolturas particulares para las rutinas como DSYRK (k actualización rango simétrico) , lo que permitiría calcular una matriz de Gram un poco más rápido que dot (MT, M).

(Parecería insignificante preocuparse por la optimización de un factor constante de 2 veces en tiempo y/o espacio, pero puede marcar la diferencia a ese umbral de la magnitud del problema que puede manejar en una sola máquina ...)

Answer 3

Ésta es pitón sencillo y no numpy, pero me acaba de lanzar juntos una rutina para llenar una matriz simétrica (y un programa de prueba para asegurarse de que es correcta):

import random 

# fill a symmetric matrix with costs (i.e. m[x][y] == m[y][x] 
# For demonstration purposes, this routine connect each node to all the others 
# Since a matrix stores the costs, numbers are used to represent the nodes 
# so the row and column indices can represent nodes 

def fillCostMatrix(dim):  # square array of arrays 
    # Create zero matrix 
    new_square = [[0 for row in range(dim)] for col in range(dim)] 
    # fill in main diagonal 
    for v in range(0,dim): 
     new_square[v][v] = random.randrange(1,10) 

    # fill upper and lower triangles symmetrically by replicating diagonally 
    for v in range(1,dim): 
     iterations = dim - v 
     x = v 
     y = 0 
     while iterations > 0: 
      new_square[x][y] = new_square[y][x] = random.randrange(1,10) 
      x += 1 
      y += 1 
      iterations -= 1 
    return new_square 

# sanity test 
def test_symmetry(square): 
    dim = len(square[0]) 
    isSymmetric = '' 
    for x in range(0, dim): 
     for y in range(0, dim): 
      if square[x][y] != square[y][x]: 
       isSymmetric = 'NOT' 
    print "Matrix is", isSymmetric, "symmetric" 

def showSquare(square): 
    # Print out square matrix 
    columnHeader = ' ' 
    for i in range(len(square)): 
     columnHeader += ' ' + str(i) 
    print columnHeader 

    i = 0; 
    for col in square: 
     print i, col # print row number and data 
     i += 1 

def myMain(argv): 
    if len(argv) == 1: 
     nodeCount = 6 
    else: 
     try: 
      nodeCount = int(argv[1]) 
     except: 
      print "argument must be numeric" 
      quit() 

    # keep nodeCount <= 9 to keep the cost matrix pretty 
    costMatrix = fillCostMatrix(nodeCount) 
    print "Cost Matrix" 
    showSquare(costMatrix) 
    test_symmetry(costMatrix) # sanity test 
if __name__ == "__main__": 
    import sys 
    myMain(sys.argv) 

# vim:tabstop=8:shiftwidth=4:expandtab 

Answer 4

Hay una serie de bien - formas conocidas de almacenar matrices simétricas para que no necesiten ocupar n^2 elementos de almacenamiento. Además, es posible reescribir las operaciones comunes para acceder a estos medios de almacenamiento revisados.El trabajo definitivo es Golub y Van Loan, Matrix Computations, 3ª edición 1996, Johns Hopkins University Press, secciones 1.27-1.2.9. Por ejemplo, al citarlos desde el formulario (1.2.2), en una matriz simétrica solo es necesario almacenar A = [a_{i,j} ] para i >= j. Entonces, suponiendo que la vector que sostiene la matriz se denota V, y que A es n-por-n, poner a_{i,j} en

V[(j-1)n - j(j-1)/2 + i] 

Esto supone 1-indexación.

Golub y Van Loan ofrecen un algoritmo 1.2.3 que muestra cómo acceder a una V almacenada para calcular y = V x + y.

Golub y Van Loan también proporcionan una forma de almacenar una matriz en forma diagonalmente dominante. Esto no ahorra almacenamiento, pero admite el acceso fácil para ciertos tipos de operaciones.

Answer 5

Es trivial completar pitónicamente [i][j] si se completa [j][i]. La cuestión del almacenamiento es un poco más interesante. Se puede aumentar la clase numpy array con un atributo packed que es útil tanto para guardar el almacenamiento como para leer posteriormente los datos.

class Sym(np.ndarray): 

    # wrapper class for numpy array for symmetric matrices. New attribute can pack matrix to optimize storage. 
    # Usage: 
    # If you have a symmetric matrix A as a shape (n,n) numpy ndarray, Sym(A).packed is a shape (n(n+1)/2,) numpy array 
    # that is a packed version of A. To convert it back, just wrap the flat list in Sym(). Note that Sym(Sym(A).packed) 


    def __new__(cls, input_array): 
     obj = np.asarray(input_array).view(cls) 

     if len(obj.shape) == 1: 
      l = obj.copy() 
      p = obj.copy() 
      m = int((np.sqrt(8 * len(obj) + 1) - 1)/2) 
      sqrt_m = np.sqrt(m) 

      if np.isclose(sqrt_m, np.round(sqrt_m)): 
       A = np.zeros((m, m)) 
       for i in range(m): 
        A[i, i:] = l[:(m-i)] 
        A[i:, i] = l[:(m-i)] 
        l = l[(m-i):] 
       obj = np.asarray(A).view(cls) 
       obj.packed = p 

      else: 
       raise ValueError('One dimensional input length must be a triangular number.') 

     elif len(obj.shape) == 2: 
      if obj.shape[0] != obj.shape[1]: 
       raise ValueError('Two dimensional input must be a square matrix.') 
      packed_out = [] 
      for i in range(obj.shape[0]): 
       packed_out.append(obj[i, i:]) 
      obj.packed = np.concatenate(packed_out) 

     else: 
      raise ValueError('Input array must be 1 or 2 dimensional.') 

     return obj 

    def __array_finalize__(self, obj): 
     if obj is None: return 
     self.packed = getattr(obj, 'packed', None) 

`` `