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:

It will have 2 fields for its eigendecomposition: U, D.
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?

Before going the whole subclass route, see if you can write several functions that do most of what you need. — hpaulj
– hpaulj, Commented Feb 3, 2024 at 16:42
@hpaulj, I can. But I want the decomposition to be cached. The examples I found go different direction from my needs. — David
– David, Commented Feb 3, 2024 at 16:57
Another option is a whole new object, with the symmetric array and its decomposition as attributes (or properties). How much of the whole array mechanisms do you need? For example I can't imagine doing much indexing and maintaining symmetry and/or the distinction between a view and copy. Or doing broadcasting, or adding or removing dimensions. If the symmetric array must always be 2d, the np.matrix subclass may be a model. np.ma is another subclass that maintains data and mask attributes while alk of the normal array things. — hpaulj
– hpaulj, Commented Feb 3, 2024 at 19:16
@hpaulj, I thought about it. Yet my in the long run I'd like to add more overloaded operations hence I think subclassing nd.array is better. Yet couldn't find any good example for adding a property. — David
– David, Commented Feb 5, 2024 at 6:28

Bob · Accepted Answer · 2024-02-06 06:57:34Z

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.

Collectives™ on Stack Overflow

Subclassing numpy array for symmetric matrices

1 Answer 1

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Comments

Your Answer

Sign up or log in

Post as a guest

Related