Numpy的“智能”对称矩阵

Question

在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

如果能有一个自动的厄米矩阵就更好了，虽然我在写这段话的时候并不需要它。

Answer 1

有很多大家都知道的方法可以存储对称矩阵，这样就不需要占用 n^2 个存储空间。而且，我们可以重新编写一些常见的操作，以便更好地使用这些新的存储方式。关于这方面的权威著作是 Golub 和 Van Loan 的《矩阵计算》，第三版，1996 年出版，约翰霍普金斯大学出版社，具体内容在 1.27-1.2.9 节。例如，他们在形式 (1.2.2) 中提到，对于一个对称矩阵，只需要存储 A = [a_{i,j} ] 当 i >= j 时。然后，假设用一个向量来表示这个矩阵，记作 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。

此外，他们还提供了一种将矩阵存储为对角主导形式的方法。这种方法虽然不节省存储空间，但可以方便某些其他类型的操作。

Answer 2

我也一直在想如何更好地处理numpy中的对称矩阵这个问题。

经过研究，我觉得答案可能是numpy在处理对称矩阵时受到了一些底层BLAS库的内存布局的限制。

虽然一些BLAS库确实利用了对称性来加快对称矩阵的计算，但它们仍然使用和完整矩阵一样的内存结构，也就是说，它们需要的空间是n^2，而不是n(n+1)/2。它们只是被告知这个矩阵是对称的，所以只使用上三角或下三角的值。

一些scipy.linalg的函数确实接受一些参数（比如在linalg.solve中使用sym_pos=True），这些参数会传递给BLAS库，虽然如果numpy能提供更多这样的支持就更好了，特别是像DSYRK（对称秩k更新）这样的函数包装，这样可以比用dot(M.T, M)计算Gram矩阵快很多。

虽然担心在时间和空间上优化一个2倍的常数因子可能显得有些挑剔，但这确实会影响到你在一台机器上能处理多大规模的问题...

Answer 3

如果你能在计算之前把矩阵对称化，下面的方法应该会比较快：

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…