Numpy“聪明”的对称matrix
在[i][j]被写入时，numpy中是否有一个智能且空间对称的matrix自动（并且透明地）填充[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 == aT) # for any symmetric matrix
一个自动的埃米特文也会很好，虽然在写作的时候我不需要这个。
如果你能够在计算之前对matrix进行对称化，那么下面的内容应该是相当快的：
def symmetrize(a): return a + aT - numpy.diag(a.diagonal())
这在合理的假设下工作（比如在运行symmetrize之前，不要同时a[0, 1] = 42和矛盾a[1, 0] = 123 ）。
如果你真的需要一个透明的对称化，你可以考虑子类化numpy.ndarray，并简单地重新定义__setitem__ ：
class SymNDArray(numpy.ndarray): def __setitem__(self, (i, j), value): super(SymNDArray, self).__setitem__((i, j), value) super(SymNDArray, self).__setitem__((j, i), value) def symarray(input_array): """ Returns a symmetrized version of the array-like input_array. Further assignments to the array are 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!
（或者等同于matrix而不是数组，取决于您的需要）。 这种方法甚至可以处理更复杂的赋值，如a[:, 1] = -1 ，它正确设置a[1, :]元素。
请注意，Python 3删除了编写def …(…, (i, j),…)的可能性，因此在使用Python 3运行之前，必须稍微调整代码： def __setitem__(self, indexes, value): (i, j) = indexes …
numpy中对称matrix最优化处理的一般问题也给我带来了困扰。
看了之后，我认为答案可能是numpy有限的内存布局支持BLAS对称matrix的基础程序。
尽pipe一些BLAS例程利用对称性来加速对称matrix上的计算，但它们仍然使用与全matrix相同的存储器结构，即n^2空间而不是n(n+1)/2 。 只是他们被告知matrix是对称的，只能使用上三angular或下三angular的值。
一些scipy.linalg例程接受标志（比如sym_pos=True上的sym_pos=True ），这些标志被传递给BLAS例程，尽pipe在numpy中支持这种标志会更好，特别是DSYRK（对称rank k更新），这将允许计算格点matrix比点（MT，M）快得多。
（可能看起来太麻烦了，担心在时间和/或空间上优化2倍的常量因素，但是它可以在单个机器上能够pipe理的问题的阈值上有所改变…）
有许多众所周知的存储对称matrix的方法，所以它们不需要占用n ^ 2个存储元件。 而且，重写通用操作来访问这些修改后的存储方法是可行的。 最后的工作是Golub和Van Loan， Matrix Computations ，1996年第3版，约翰霍普金斯大学出版社，第1.27-1.2.9节。 例如，从forms（1.2.2）中引用它们，在对称matrix中只需要存储A = [a_{i,j} ] ， i >= j 。 那么，假设保持matrix的向量被表示为V，并且A是n乘n，则将a_{i,j}放入
V[(j-1)n - j(j-1)/2 + i]
这假设1索引。
Golub和Van Loan提供了一个algorithm1.2.3，它演示了如何访问这样存储的V来计算y = V x + y 。
Golub和Van Loan也提供了一种以对angular占优forms存储matrix的方法。 这不会节省存储空间，但支持某些其他types的操作。
这是普通的python，而不是numpy，但我只是把一个例程，以填补对称matrix（和testing程序，以确保它是正确的）：
import random # fill a symmetric matrix with costs (ie 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
如果[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[:(mi)] A[i:, i] = l[:(mi)] l = l[(mi):] 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)
“`
