制作一个矩阵对称，就地与不在场

Question

我有一个矩阵，它应该是对称的（它是对称的逆），但它不完全是由于反演中的数字错误等。因此，我添加了一个使矩阵对称的步骤（由a = .5(a+a'开始），如果我在原地进行（不合适的地方），我会看到一个数字灾难。代码：

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：

换位是一个视图。

作为一个结果，下面的代码是错误的，不会让一个矩阵对称：

a += a.T 

它将为小数组（因为缓冲）工作，但失败的大一点，在不可预知的方式。

的原因是，在同一时间a正在用a.T更新，a.T实际改变（因为它是一个memoryview的a），并从而更新的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的临时副本。