CODEWITHSUNDEEP
pythonarraysnumpyscipy
blaz
blaz

Reputation: 4248

Inverse of a symmetric matrix

An inverse of a real symmetric matrix should in theory return a real symmetric matrix (the same is valid for Hermitian matrices). However, when I compute the inverse with numpy or scipy the returned matrix is asymmetric. I understand that this is due to numerical error.

What is the best way to avoid this asymmetry? I want it to be valid mathematically, in order that it does not propagate the error further when I use it in my computations.

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]

Upvotes: 3

Views: 2969

Answers (2)

user3684792
user3684792

Reputation: 2611

Luckily for you, this inverse is symmetric. Unluckily for you you can't compare floating points this way:

>>> 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

Luckily for you, you can use numpy all close 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 its your lucky day

Edit: Wow, I am quite surprised cels answer looks to be 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]

Upvotes: 1

cel
cel

Reputation: 31349

This simple change should convince you, that the inverse is indeed a symmetric matrix. Not mathematically, but at least numerically - that is up to a small error threshold 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)))

Which outputs:

Inverse of matrix A is symmetric

Upvotes: 1

Related Questions