python - 对称矩阵的逆

Question

理论上，实数对称矩阵的逆应该返回实数对称矩阵(这同样适用于 Hermitian 矩阵)。但是，当我用 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

这个简单的改变应该让您相信逆矩阵确实是一个对称矩阵。不是数学上的，但至少是数字上的——这是一个小的错误阈值 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)))

哪些输出:

Inverse of matrix A is symmetric