Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Numpy ‘smart’ symmetric matrix

Is there a smart and space-efficient symmetric matrix in numpy which automatically (and transparently) fills the position at [j][i] when [i][j] is 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 == a.T) # for any symmetric matrix

An automatic Hermitian would also be nice, although I won’t need that at the time of writing.

like image 916
Debilski Avatar asked Apr 03 '10 22:04

Debilski


People also ask

How do you check if a matrix is symmetric in Numpy?

A Simple solution is to do following. 1) Create transpose of given matrix. 2) Check if transpose and given matrices are same or not.

How do you define a symmetric matrix in python?

A matrix is called symmetric if a_{ij} is equal to a_{ji}. Thanks to this rule, an N \times N symmetric matrix needs to store only (N + 1) \cdot \frac{N}{2} elements instead of N^2 elements needed to be stored in case of a classic matrix.

How do you check if a matrix is symmetric or not in Python?

A square matrix is said to be symmetric matrix if the transpose of the matrix is same as the given matrix. Symmetric matrix can be obtain by changing row to column and column to row. Recommended: Please try your approach on {IDE} first, before moving on to the solution.

How do you know if a matrix is symmetric?

How to check Whether a Matrix is Symmetric or Not? Step 1- Find the transpose of the matrix. Step 2- Check if the transpose of the matrix is equal to the original matrix. Step 3- If the transpose matrix and the original matrix are equal, then the matrix is symmetric.


2 Answers

If you can afford to symmetrize the matrix just before doing calculations, the following should be reasonably fast:

def symmetrize(a):     """     Return a symmetrized version of NumPy array a.      Values 0 are replaced by the array value at the symmetric     position (with respect to the diagonal), i.e. if a_ij = 0,     then the returned array a' is such that a'_ij = a_ji.      Diagonal values are left untouched.      a -- square NumPy array, such that a_ij = 0 or a_ji = 0,      for i != j.     """     return a + a.T - numpy.diag(a.diagonal()) 

This works under reasonable assumptions (such as not doing both a[0, 1] = 42 and the contradictory a[1, 0] = 123 before running symmetrize).

If you really need a transparent symmetrization, you might consider subclassing numpy.ndarray and simply redefining __setitem__:

class SymNDArray(numpy.ndarray):     """     NumPy array subclass for symmetric matrices.      A SymNDArray arr is such that doing arr[i,j] = value     automatically does arr[j,i] = value, so that array     updates remain symmetrical.     """      def __setitem__(self, (i, j), value):         super(SymNDArray, self).__setitem__((i, j), value)                             super(SymNDArray, self).__setitem__((j, i), value)                      def symarray(input_array):     """     Return a symmetrized version of the array-like input_array.      The returned array has class SymNDArray. Further assignments to the array     are thus 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 the equivalent with matrices instead of arrays, depending on your needs). This approach even handles more complicated assignments, like a[:, 1] = -1, which correctly sets a[1, :] elements.

Note that Python 3 removed the possibility of writing def …(…, (i, j),…), so the code has to be slightly adapted before running with Python 3: def __setitem__(self, indexes, value): (i, j) = indexes

like image 189
Eric O Lebigot Avatar answered Sep 30 '22 07:09

Eric O Lebigot


The more general issue of optimal treatment of symmetric matrices in numpy bugged me too.

After looking into it, I think the answer is probably that numpy is somewhat constrained by the memory layout supportd by the underlying BLAS routines for symmetric matrices.

While some BLAS routines do exploit symmetry to speed up computations on symmetric matrices, they still use the same memory structure as a full matrix, that is, n^2 space rather than n(n+1)/2. Just they get told that the matrix is symmetric and to use only the values in either the upper or the lower triangle.

Some of the scipy.linalg routines do accept flags (like sym_pos=True on linalg.solve) which get passed on to BLAS routines, although more support for this in numpy would be nice, in particular wrappers for routines like DSYRK (symmetric rank k update), which would allow a Gram matrix to be computed a fair bit quicker than dot(M.T, M).

(Might seem nitpicky to worry about optimising for a 2x constant factor on time and/or space, but it can make a difference to that threshold of how big a problem you can manage on a single machine...)

like image 35
Matt Avatar answered Sep 30 '22 07:09

Matt