マトリックスを対称で、インプレースとアウトオブプレイスにする

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）のために、問題が表示されないこと。

EDITいくつかの詳細情報は、インプレースバージョン後a-a'を見て次式で与えられます。

Answer 1

エラーがラインa += a.Tから来ています。これは、インプレース操作の既知の問題です（私は今そう述べている文書の適切な部分を見つけることができません）が、scipy lecture notesから引用：

転置図です。結果として

、次のコードが間違っているとマトリックス対称をすることはありません。

a += a.T 

それが（なぜならバッファリングの）小配列のために働くが、予測不可能な方法で、大きなもののために失敗します。

理由は、（それがmemoryview aのあるので）、従って、誤っaのいくつかの座標を更新すると同時にaがa.Tが実際に変更され、a.Tで更新されていることです。あなたはインプレース行列を対称にしたい場合は、次の操作を行うことができ

：

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

それとも、あなたの表記法に固執する場合：copy以来

a += a.T.copy() 

が作成されます。更新されないa.Tの一時的なコピーです。