from - python matrices tutorial

Numpy matriz simétrica "inteligente"

Question

¿Hay una matriz simétrica inteligente y eficiente en el espacio en numpy que automáticamente (y transparentemente) 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 hermitiano automático también sería agradable, aunque no necesitaré eso al momento de escribirlo.

Answer 1

El problema más general del tratamiento óptimo de las matrices simétricas en numpy también me molestó.

Después de analizarlo, creo que la respuesta probablemente sea que numpy está algo restringido por el diseño de memoria soportado 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.

Algunas de las rutinas scipy.linalg aceptan banderas (como sym_pos=True en linalg.solve ) que pasan a las rutinas BLAS, aunque sería bueno contar con más soporte para esto en numpy, en particular envoltorios para rutinas como DSYRK (rango simétrico k actualización), que permitiría calcular una matriz de Gram un poco más rápido que el punto (MT, M).

(Puede parecer quisquilloso preocuparse por la optimización de un factor constante de 2 veces en tiempo y / o espacio, pero puede marcar la diferencia en ese umbral de cuán grande es el problema que puede manejar en una sola máquina ...)

Answer 2

Es trivial completar pitónicamente [i][j] si [j][i] está lleno. La pregunta de almacenamiento es un poco más interesante. Se puede aumentar la clase de matriz numpy con un atributo packed que es útil tanto para guardar el almacenamiento como para leer los datos más tarde.

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)

`` `

Answer 3

Esto es simple y simple, pero simplemente armé una rutina para llenar una matriz simétrica (y un programa de prueba para asegurarme de que sea 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 formas bien 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 de la forma (1.2.2), en una matriz simétrica solo se necesita almacenar A = [a_{i,j} ] para i >= j . Luego, suponiendo que el vector que contiene la matriz se denota V, y que A es n-by-n, ponga a 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

Si puede permitirse simétrizar la matriz justo antes de hacer los cálculos, lo siguiente debe ser razonablemente rápido:

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

Esto funciona bajo suposiciones razonables (como no hacer tanto a[0, 1] = 42 y la contradictoria a[1, 0] = 123 antes de ejecutar symmetrize ).

Si realmente necesita una simetría transparente, puede considerar la creación de subclases de numpy.ndarray y simplemente redefinir __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, según sus necesidades). Este enfoque incluso maneja asignaciones más complicadas, como a[:, 1] = -1 , que establece correctamente a[1, :] .

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