在 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 也不错，尽管在撰写本文时我不需要它。
【问题讨论】：
如果您有能力在进行计算之前使矩阵对称，那么以下应该相当快：
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...
【讨论】：
print 时，将其写为等效的__getitem__(self, (i, j)) 会失败。原因是print 使用整数索引调用__getitem__()，因此即使是简单的print 也需要做更多的工作。 __setitem__() 的解决方案适用于print（显然），但遇到类似的问题：a[0] = [1, 2, 3] 不起作用，出于同样的原因（这不是一个完美的解决方案）。 __setitem__() 解决方案的优点是更健壮，因为内存中的数组是正确的。还不错。 :)
在 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 倍常数因子似乎有些挑剔，但它可能会影响您在单台机器上可以处理多大问题的阈值...）
【讨论】：
有许多众所周知的存储对称矩阵的方法，因此它们不需要占用 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 还提供了一种以对角占优形式存储矩阵的方法。这不会节省存储空间，但支持某些其他类型的操作的就绪访问。
【讨论】：
这是普通的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
【讨论】：
要构造一个沿主对角线对称且主对角线上为 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]))
这是一种特殊情况，但最近我使用这种矩阵来表示网络邻接关系。
希望对您有所帮助。 干杯。
【讨论】：
如果填写[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)
```
【讨论】：