对对称矩阵的 numpy 数组进行子类化

Question

我想创建一个专门用于对称矩阵的 numpy 子类。 我希望它就像任何 numpy 数组一样，只有 2 个差异：

1. 它将有 2 个用于特征分解的字段：

   U

   、

   D

   。
2. 在施工时它将强制对称。
   就像

   my_symmetric_arr() = 0.5*(m+m.T)

   一样，其中

   m

   是由默认构造函数创建的常规矩阵。

怎样才能以最简单的方式做到这一点？

Answer 1

对 numpy 数组进行子类化可能很棘手，但幸运的是您只是偶尔做一次。 Numpy 文档提供了有关该主题的良好资源there

这里我给你一个例子，你可以从输入数组构建对称数组，当你分配给数组的切片时，它也会分配转置的切片

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)
        

这里是一个用法示例

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)

打印以下数组

[[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.]]

对于特征分解，我不会将其保留为数组的属性，numpy 已经为对称特征分解提供了 eigh。