对称矩阵的逆矩阵

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]

- blaz

---

问题出在这一行：if (a_symm_inv == a_symm_inv.T).all():。浮点数的算术运算并不是精确的，因此您不应该尝试进行精确的比较。 - cel

但我不明白你的问题。既然你知道结果是对称矩阵，并且想要强制浮点数相等（这可能没有太多意义），你可以用上三角矩阵的值覆盖下三角矩阵。 - cel

我希望它在数学上是有效的。你不能用真正的计算机计算出精确的逆矩阵。有几种数值算法可以近似求解逆矩阵。我认为你想得太多了。误差源于计算机浮点运算的限制。在一个完美的世界里，你拥有无限精度的处理器，你会得到一个完美的结果 :) - cel

@blaz：也许你想让矩阵保持对称性有一个非常具体的原因，但如果你关注的是数值精度，那么你并不能通过这种方式获得更精确的结果。你只是更喜欢上三角矩阵的误差。据我所知，你的请求没有意义。 - cd98

我同意所有的回复。从数值角度来看，这个错误是可以预料的。正如@cd98所提到的，我更喜欢上三角矩阵的误差。我只是在想，也许numpy中已经实现了识别矩阵对称性的功能，因此只计算逆矩阵的下（或上）三角形，并将其复制到另一个三角形中。在我看来，这样会更快。 - blaz

Answer 1

幸运的是，这个倒数是对称的。不幸的是，你不能用这种方式比较浮点数：

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

幸运的是，你可以使用numpy all close来解决这个问题。http://docs.scipy.org/doc/numpy/reference/generated/numpy.allclose.html

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

看起来今天是你的幸运日

编辑：哇，我很惊讶Cels的答案似乎比这个更快：

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

- user3684792

---

我认为对称矩阵的逆矩阵总是对称的。 - Kyle_S-C

Answer 2

这个简单的变化应该让你相信，逆矩阵确实是一个对称矩阵。虽然不是数学上的，但至少在数值上 - 即在一个小误差阈值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

- cel

---

使用numpy all close - 它比手动实现更好，而且几乎完全相同。 - user3684792

@user3684792 我不明白为什么这样会“更好”。它掩盖了我们正在使用一个 epsilon 的事实。 - cel

你说得对，我对它们进行了分析，结果发现numpy慢了4倍。看看吧。 - user3684792