Numpy“智能”对称矩阵

Posted 2023-02-23

【中文标题】Numpy“智能”对称矩阵

【英文标题】：Numpy ‘smart’ symmetric matrix

【发布时间】：2011-02-04 02:24:39

【问题描述】：

在 numpy 中是否有一个智能且节省空间的对称矩阵，当 [i][j] 被写入时，它会自动（并且透明地）填充 [j][i] 的位置？

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

自动 Hermitian 也不错，尽管在撰写本文时我不需要它。

【问题讨论】：

如果可以解决您的问题，您可以考虑将答案标记为已接受。 :)

我想等待一个更好的（即内置和内存效率高的）答案的出现。当然，你的回答没有问题，所以我还是会接受。

【参考方案1】：

如果您有能力在进行计算之前使矩阵对称，那么以下应该相当快：

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())

这在合理的假设下有效（例如在运行 symmetrize 之前不同时执行 a[0, 1] = 42 和矛盾的 a[1, 0] = 123）。

如果您真的需要透明的对称化，您可以考虑继承 numpy.ndarray 并简单地重新定义 __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!

（或使用矩阵代替数组的等价物，具体取决于您的需要）。这种方法甚至可以处理更复杂的赋值，例如 a[:, 1] = -1，它正确设置了 a[1, :] 元素。

请注意，Python 3 消除了编写 def …(…, (i, j),…) 的可能性，因此在使用 Python 3 运行之前必须稍微修改代码：def __setitem__(self, indexes, value): (i, j) = indexes...

【讨论】：

实际上，如果你对它进行子类化，你不应该覆盖 setitem，而应该覆盖 getitem，这样你就不会在创建矩阵。

这是一个非常有趣的想法，但是当在子类实例数组上执行简单的print 时，将其写为等效的__getitem__(self, (i, j)) 会失败。原因是print 使用整数索引调用__getitem__()，因此即使是简单的print 也需要做更多的工作。 __setitem__() 的解决方案适用于print（显然），但遇到类似的问题：a[0] = [1, 2, 3] 不起作用，出于同样的原因（这不是一个完美的解决方案）。 __setitem__() 解决方案的优点是更健壮，因为内存中的数组是正确的。还不错。 :)

你的建议听起来像blog.sopticek.net/2016/07/24/…...你确认它几乎一样吗？麻烦的是这优化了内存使用，而不是计算时间。我正在寻找 python 方法来加速对称矩阵上的一些简单计算。如果您有信息，请告诉我。

这个答案不会节省内存，因此与引用链接中的方法有很大不同。现在，使用对称矩阵节省时间通常涉及使用专用算法而不是通用算法，例如在 NumPy 中使用 eigh() 而不是 eig()。

【参考方案2】：

在 numpy 中优化对称矩阵的更普遍问题也困扰着我。

在研究之后，我认为答案可能是 numpy 在某种程度上受到对称矩阵的底层 BLAS 例程支持的内存布局的限制。

虽然一些 BLAS 例程确实利用对称性来加速对称矩阵的计算，但它们仍然使用与完整矩阵相同的内存结构，即 n^2 空间而不是 n(n+1)/2。只是他们被告知矩阵是对称的，并且只能使用上三角形或下三角形中的值。

一些scipy.linalg 例程确实接受标志（如linalg.solve 上的sym_pos=True），这些标志会传递给BLAS 例程，尽管在numpy 中对此提供更多支持会很好，特别是像DSYRK 这样的例程的包装器（对称秩 k 更新），这将允许 Gram 矩阵的计算比 dot(MT, M) 快一点。

（担心在时间和/或空间上优化 2 倍常数因子似乎有些挑剔，但它可能会影响您在单台机器上可以处理多大问题的阈值...）

【讨论】：

问题是关于如何通过分配单个条目来自动创建对称矩阵（而不是关于如何指示 BLAS 在其计算中使用对称矩阵或原则上如何存储对称矩阵更有效）。

这个问题也是关于空间效率的，所以 BLAS 问题是热门话题。

@EOL，问题不在于如何通过单个条目的分配自动创建对称矩阵。

当然，“创建”可以更恰当地替换为“更新”。现在，由于问题明确是关于在设置 M_ji 时透明地设置 M_ji ，而这个答案与此无关，你明白这本质上就是我提出的观点。问题是关于如何有效地做到这一点（而不是有效地处理对称矩阵，即使这可能是正确的问题：更好地放在 cmets 上，或者作为解决更普遍问题的答案而不是仅仅讨论它)。

【参考方案3】：

有许多众所周知的存储对称矩阵的方法，因此它们不需要占用 n^2 个存储元素。此外，重写常用操作以访问这些修改后的存储方式是可行的。权威著作是 Golub 和 Van Loan，矩阵计算，1996 年第 3 版，约翰霍普金斯大学出版社，第 1.27-1.2.9 节。例如，从表格（1.2.2）中引用它们，在对称矩阵中只需要存储A = [a_i,j ] fori >= j。然后，假设保存矩阵的 vector 表示为 V，并且 A 是 n×n，将 a_i,j 放入

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

这假定为 1 索引。

Golub 和 Van Loan 提供了一个算法 1.2.3，它展示了如何访问这样一个存储的 V 来计算 y = V x + y。

Golub 和 Van Loan 还提供了一种以对角占优形式存储矩阵的方法。这不会节省存储空间，但支持某些其他类型的操作的就绪访问。

【讨论】：

还有 Rectangular Full Packed storage (RFP)，例如 Lapack ZPPTRF 使用它。 numpy 支持吗？

@isti_spl：不，但你可以实现一个包装器

【参考方案4】：

这是普通的python而不是numpy，但我只是拼凑了一个例程来填充 一个对称矩阵（以及一个确保它正确的测试程序）：

import random

# fill a symmetric matrix with costs (i.e. 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

【参考方案5】：

要构造一个沿主对角线对称且主对角线上为 0 的 NxN 矩阵，您可以这样做：

a = np.array([1, 2, 3, 4, 5])
b = np.zeros(shape=(a.shape[0], a.shape[0]))
upper = np.triu(b + a)
lower = np.tril(np.transpose(b + a))
D = (upper + lower) * (np.full(a.shape[0], fill_value=1) - np.eye(a.shape[0]))

这是一种特殊情况，但最近我使用这种矩阵来表示网络邻接关系。

希望对您有所帮助。 干杯。

【参考方案6】：

如果填写[j][i]，用Python 方式填写[i][j] 很简单。存储问题有点意思。可以使用packed 属性来扩充 numpy 数组类，该属性对于节省存储空间和以后读取数据都很有用。

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

```