python - 使矩阵对称、就地与异地

Question

我有一个应该是对称的矩阵(它是对称的逆矩阵)，但它并不完全是由于反演等中的数值误差造成的。

因此，我添加了使矩阵对称的步骤(通过 a = .5(a+a'))，如果我就地执行此操作(out- of-place 就可以了)。代码:

import numpy as np

def check_sym(x):
    print("||a-a'||^2 = %e" % np.sum((x - x.T)**2))

# make a symmetric matrix
dim = 100
a = np.random.randn(dim,dim)
a = np.matmul(a, a.T)
b = a.copy()

check_sym(a)

print("symmetrizing in-place")
a += a.T
a *= .5
check_sym(a)

print("symmetrizing out-of-place")
b = .5 * (b + b.T)
check_sym(b)

输出是:

||a-a'||^2 = 1.184044e-26
symmetrizing in-place
||a-a'||^2 = 7.313593e+04
symmetrizing out-of-place
||a-a'||^2 = 0.000000e+00

请注意，对于较低维度(例如 dim=10)，该问题不会出现。

编辑通过查看就地版本后的a-a'给出了更多信息:

Answer 1

错误来自行a += a.T。这是就地操作的一个已知问题(我现在找不到正确的文档来说明这一点)，但引用自 scipy lecture notes :

The transposition is a view.

As a results, the following code is wrong and will not make a matrix symmetric:

a += a.T

It will work for small arrays (because of buffering) but fail for large one, in unpredictable ways.

原因是，在使用 a.T 更新 a 的同时，a.T 实际上也在发生变化(因为它是一个 a 的 >memoryview)，从而错误地更新了 a 的某些坐标。

如果您想就地对称化矩阵，您可以执行以下操作:

a = np.random.rand(4,4)
a[np.tril_indices_from(a)] = a.T[np.tril_indices_from(a)]

或者，如果您想坚持您的符号:

a += a.T.copy()

因为 copy 将创建 a.T 的临时副本，该副本不会更新。