Numpy的 "智能 "对称矩阵

Question 1

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也会很好，尽管在写这篇文章时我不需要这个。

Question 2


          
           如果你能在进行计算前对矩阵进行对称，下面的计算应该是相当快的。
          
          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...

Question 3


          
           
            
             在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)快得多。T, M）。
            
            
             (担心对时间和/或空间的2倍恒定系数进行优化似乎有些吹毛求疵，但这对你在单台机器上能处理多大的问题的门槛是有影响的。)

Question 4


          
           
            
             有一些众所周知的存储对称矩阵的方法，因此它们不需要占用n^2个存储元素。  此外，重写常见的操作来访问这些修改后的存储方式也是可行的。  最权威的工作是Golub和Van Loan。
             
              矩阵计算
             
             ，1996年第三版，约翰霍普金斯大学出版社，第1.27-1.2.9节。  例如，从表格（1.2.2）中引用它们，在一个对称矩阵中只需要存储
             
              A = [a_{i,j} ]
             
             ，用于
             
              i >= 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还提供了一种以对角线主导形式存储矩阵的方法。  这并不节省存储空间，但支持对某些其他类型操作的随时访问。

Question 5


          
           
            
             
              这是普通的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:
            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

Question 6


          
           
            
             
              
               为了构建一个沿主对角线对称的NxN矩阵，并且在主对角线上有0，你可以做:
              
              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]))
这是一种特殊情况，但最近我把这种矩阵用于网络邻接表示。
希望这有帮助。

Question 7


          
           
            
             
              
               
                如果
                
                 [j][i]
                
                被填上了，那么用Python填上
                
                 [i][j]
                
                是很容易的。  存储问题更有趣一些。我们可以在numpy数组类中增加一个
                
                 packed
                
                的属性，这对于节省存储空间和以后读取数据都很有用。
               
               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):