Subclassing numpy array for symmetric matrices

Question

I want to create a numpy subclass specialized for symmetric matrices. I want it to be just like any numpy array with only 2 differences:

1. It will have 2 fields for its eigendecomposition: U, D.
2. On construction it will force symmetry.
   Like my_symmetric_arr() = 0.5*(m+m.T) where m is the regular matrix created by the default constructor.

How can that be done in the simplest way?

Answer 1

Subclassing numpy arrays can be tricky, but fortunately is something you do only once in a while. Numpy documentation provides a good resource on the topic there

Here I give you an example from which you can build a symmetric array from an input array, and when you assign to slices of the array it will also assign the transposed slice

import numpy as np

class SymmetricArray(np.ndarray):
    def __new__(cls, input_array):
        '''
        Creates SymmetricArray from the symmetric part of the input array
        (input_array + input_array.T) / 2
        If the array has more than two dimensions the two last axes are
        treated as matrix dimensions, and the rest will be broadcast indices
        '''
        assert input_array.shape[-1] == input_array.shape[-2]
        transposed_input_array = input_array.transpose(list(range(input_array.ndim - 2)) + [-1, -2])
        return (0.5 * (input_array + transposed_input_array)).view(SymmetricArray)

    def __setitem__(self, loc, value):
        '''
        Assigns to slices or items in the array by preserving symmetry
        of the last two dimensions
        '''
        if type(loc) != tuple:
            loc = (loc,)
        if len(loc) < self.ndim:
            loc += (self.ndim - len(loc)) * (slice(None, None, None),)
        
        loc_t = loc[:-2] + (loc[-1], loc[-2])
        value = np.asarray(value)
        value_t = np.asarray(value)
        if value.ndim >= 2:
            value_t = value.transpose(
                       list(range(value.ndim - 2)) + [-1, -2])
        super().__setitem__(loc, value)
        super().__setitem__(loc_t, value_t)
        

Here an example usage

rng = np.random.default_rng(0)
a = SymmetricArray(rng.integers(0, 50, (4,4))*2)
print(a)
a[:, 1] = 0
print(a)
a[3,1] = -4
print(a)
a[2:, :2] = 11
print(a)

That prints the following arrays

[[84. 46. 33. 38.]
 [46.  4. 43. 30.]
 [33. 43. 64. 93.]
 [38. 30. 93. 72.]]

[[84.  0. 33. 38.]
 [ 0.  0.  0.  0.]
 [33.  0. 64. 93.]
 [38.  0. 93. 72.]]

[[84.  0. 33. 38.]
 [ 0.  0.  0. -4.]
 [33.  0. 64. 93.]
 [38. -4. 93. 72.]]

[[84.  0. 11. 11.]
 [ 0.  0. 11. 11.]
 [11. 11. 64. 93.]
 [11. 11. 93. 72.]]

For the eigendecomposition I wouldn't keep that as an attribute of the array, numpy already provides the eigh for symmetric eigendecomposition.