Inverse symmetric matrix
The inverse of a real symmetric matrix should theoretically return a real symmetric matrix (the same is true for Hermitian matrices). However, when I calculate the inverse with numpy or scipy, the returned matrix is asymmetric. I understand this is due to a numeric error.
What is the best way to avoid this asymmetry? I want it to be valid mathematically so that it doesn't propagate the error further when I use it in my calculations.
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)))
EDIT
This is my solution to the problem:
math_symm = (np.triu_indices(len(a_symm_inv), 1))
a_symm_inv[math_symm]=np.tril(a_symm_inv, -1).T[math_symm]
source to share
Luckily for you, this reverse is symmetrical. Unfortunately for you, you cannot compare floats like this:
>>> 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
Lucky for you, you can use numpy to solve this problem http://docs.scipy.org/doc/numpy/reference/generated/numpy.allclose.html
>>> np.allclose(a_symm_inv, a_symm_inv_T)
True
Looks like your lucky day
Edit: Wow, I'm pretty surprised that cels' answer looks faster than this:
>>> 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]
source to share
This simple change should convince you that the inverse is indeed a symmetric matrix. Not mathematically, but at least numerically - this is up to the threshold of error epsilon
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)))
What are the outputs:
Inverse of matrix A is symmetric
source to share