製作一個矩陣對稱，就地與不在場

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的臨時副本。