Smart Numpy Symmetric Matrix

Question

Is there an intelligent and spatially efficient symmetric matrix in numpy that automatically (and transparently) fills the position in [j][i] when [i][j] written to?

 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 == aT) # for any symmetric matrix 

Automatic Hermitian will also be pleasant, although I will not need it at the time of writing.

Answer 1

If you can afford to symmetry the matrix just before performing the calculations, the following should be fast enough:

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

This works under reasonable assumptions (for example, not making either a[0, 1] = 42 or the inconsistent a[1, 0] = 123 before running symmetrize ).

If you really need transparent symmetrization, you can consider subclassing numpy.ndarray and simply override __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! 

(or equivalent with matrices instead of arrays, depending on your needs). This approach even handles more complex assignments, such as a[:, 1] = -1 , which sets the elements a[1, :] a[:, 1] = -1 correctly.

Note that Python 3 removed the ability to write def …(…, (i, j),…) , so before starting with Python 3, the code

Answer 2

The more general question of optimal treatment for symmetric matrices in numpy also listened to me.

Looking at it, I think the answer is probably that numpy is somewhat limited to supporting the layout of the memory under the basic BLAS procedures for symmetric matrices.

Although some BLAS procedures use symmetry to speed up computations on symmetric matrices, they still use the same memory structure as the full matrix, i.e. n^2 space, not n(n+1)/2 . They are simply told that the matrix is ​​symmetrical and uses only the values ​​in the upper or lower triangle.

Some of the scipy.linalg routines accept flags (e.g. sym_pos=True on linalg.solve ) that are passed to BLAS procedures, although much support for this in numpy will be nice, in particular wrappers for routines such as DSYRK (symmetric rank update k ), which would allow us to calculate the Gram matrix faster than the point (MT, M).

(It may seem odd to worry about optimizing for a 2x constant factor over time and / or space, but it may affect this threshold of how big a problem you can handle on one machine ...)

Answer 3

There are a number of well-known methods for storing symmetric matrices, so they do not need to occupy n ^ 2 storage elements. You can also rewrite common operations to access these revised storage facilities. The final work is Golub and Van Loan, Matrix Computations, 3rd Edition 1996, Johns Hopkins University, sections 1.27-1.2.9. For example, quoting them from the form (1.2.2), in the symmetric matrix we need to store A = [a_{i,j} ] for i >= j . Then, assuming the matrix holding vector is denoted by V and A is n-by-n, put a_{i,j} in

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

This involves 1-indexing.

Golub and Van Loan propose Algorithm 1.2.3, which shows how to access such a stored V in order to calculate y = V x + y .

Golub and Van Loan also provide a way to store the matrix in a diagonal dominant form. This does not save memory, but maintains access to them for some other kinds of operations.

Answer 4

This is simple python, but not numpy, but I just compiled a routine to populate the symmetric matrix (and a test program to make sure it is correct):

 import random # fill a symmetric matrix with costs (ie 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 5

This is trivial for Pythonically populating [i][j] if populating [j][i] . The storage issue is a bit more interesting. You can increase the numpy array class using the packed attribute, which is useful both for storing storage and for reading data later.

 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[:(mi)] A[i:, i] = l[:(mi)] l = l[(mi):] 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) 

`` ``