대칭 행렬의 역함수

Question

실제 대칭 행렬의 역함수는 이론상 실제 대칭 행렬을 반환해야합니다 (동일한 것은 은닉 행렬에 유효합니다). 그러나 numpy 또는 scipy를 사용하여 역함수를 계산할 때 반환 된 행렬은 비대칭입니다. 이 오류는 숫자 오류로 인한 것임을 이해합니다.

이 비대칭을 피하는 가장 좋은 방법은 무엇입니까? 나는 그것을 수학적으로 유효하게하고 싶다. 계산을 위해서 사용할 때 오류를 더 전파하지 않기를 바란다.

import numpy as np 

n = 1000 
a =np.random.rand(n, n) 
a_symm = (a+a.T)/2 

a_symm_inv = np.linalg.inv(a_symm) 

if (a_symm_inv == a_symm_inv.T).all(): 
    print("Inverse of matrix A is symmetric") # This does not happen! 
else: 
    print("Inverse of matrix A is asymmetric") 
    print("Max. asymm. value: ", np.max(np.abs((a_symm_inv-a_symm_inv.T)/2))) 

편집

이 문제에 대한 내 솔루션입니다 :

math_symm = (np.triu_indices(len(a_symm_inv), 1)) 
a_symm_inv[math_symm]=np.tril(a_symm_inv, -1).T[math_symm] 

Answer 1

다행히 당신을 위해,이 역은 대칭입니다. 당신은 http://docs.scipy.org/doc/numpy/reference/generated/numpy.allclose.html

>>> np.allclose(a_symm_inv, a_symm_inv_T) 
True 

은 당신의 행운의 날

것 같은데이 문제를 해결하기 위해 모든 가까이 NumPy와 사용할 수 있습니다, 당신을 위해 다행히

>>> import numpy as np 
>>> 
>>> n = 1000 
>>> a =np.random.rand(n, n) 
>>> a_symm = (a+a.T)/2 
>>> 
>>> a_symm_inv = np.linalg.inv(a_symm) 
>>> a_symm_inv_T = a_symm_inv.T 
>>> print a_symm_inv[2,1] 
0.0505944152801 
>>> print a_symm_inv_T[2,1] 
0.0505944152801 
>>> print a_symm_inv[2,1] == a_symm_inv_T[2,1] 
False 

: 불행하게도 당신을 위해 당신이 지점이 방법을 떠 비교할 수 없습니다

편집 :

>>> import timeit 
>>> setup = """import numpy as np 
... a = np.random.rand(1000, 1000) 
... b = np.random.rand(1000, 1000) 
... def cool_comparison_function(matrix1, matrix2): 
...  epsilon = 1e-9 
...  if (np.abs(matrix1 - matrix2) < epsilon).all(): 
...    return True 
...  else: 
...    return False 
... """ 
>>> timeit.Timer("cool_comparison_function(a,b)",setup).repeat(1, 1000) 
[2.6709160804748535] 
>>> timeit.Timer("np.allclose(a,b)",setup).repeat(1, 1000) 
[11.295115947723389] 

Answer 2

: 와우, 나는 매우 놀랄 CELS 대답입니다 빨리이보다 보인다

이 간단한 변경은 역 분모가 실제로 대칭 행렬임을 확신시켜야합니다. 아니 수학적으로, 그러나 적어도 수치 - 출력하는 작은 오류 임계 값 엡실론

n = 1000 
a =np.random.rand(n, n) 
a_symm = (a+a.T)/2 

a_symm_inv = np.linalg.inv(a_symm) 
epsilon = 1e-9 
if (np.abs(a_symm_inv - a_symm_inv.T) < epsilon).all(): 
    print("Inverse of matrix A is symmetric") 
else: 
    print("Inverse of matrix A is asymmetric") 
    print("Max. asymm. value: ", np.max(np.abs((a_symm_inv-a_symm_inv.T)/2))) 

까지입니다 :

Inverse of matrix A is symmetric