如何从 numpy.dot(A, A_inv) 得出单位矩阵答案

【问题标题】：How to arrive at the unit matrix from numpy.dot(A, A_inv)如何从 numpy.dot(A, A_inv) 得出单位矩阵
【发布时间】：2013-03-20 09:20:04
【问题描述】：

我准备了一个随机数矩阵，计算它的逆矩阵并将其与原始矩阵相乘。这在理论上给出了单位矩阵。我怎样才能让numpy 为我这样做？

import numpy

A = numpy.zeros((100,100))
E = numpy.zeros((100,100))

size = 100

for i in range(size):
    for j in range(size):
        A[i][j]+=numpy.random.randint(10)
        if i == j:
            E[i][j]+=1

A_inv = numpy.linalg.linalg.inv(A)
print numpy.dot(A, A_inv)

运行代码产生

[me]machine @ numeric $ python rand_diag.py 
[[  1.00000000e+00  -7.99360578e-15  -1.14491749e-16 ...,   3.81639165e-17
   -4.42701431e-15   1.17961196e-15]
[ -5.55111512e-16   1.00000000e+00  -2.22044605e-16 ...,  -3.88578059e-16
    1.33226763e-15  -8.32667268e-16]

很明显，结果是一个单位矩阵，但并不精确，所以print numpy.dot(A, A_inv) == E 显然给出了False。我这样做是为了练习线性代数并试图找到我的机器达到其极限的矩阵大小。获得True 将具有教学吸引力。

编辑：

设置size=10000，我内存不足

[me]machine @ numeric $ Python(794) malloc:
***mmap(size=800002048) failed (error code=12)
*** error: can\'t allocate region
*** set a breakpoint in malloc_error_break to debug
Traceback (most recent call last):
File "rand_diag.py", line 14, in <module>     A_inv = numpy.linalg.linalg.inv(A)
File "/Library/Frameworks/Python.framework/Versions/7.2/lib/python2.7/site-packages/numpy/linalg/linalg.py", line 445, in inv
return wrap(solve(a, identity(a.shape[0], dtype=a.dtype)))
File "/Library/Frameworks/Python.framework/Versions/7.2/lib/python2.7/site-packages/numpy/linalg/linalg.py", line 323, in solve
a, b = _fastCopyAndTranspose(t, a, b)
File "/Library/Frameworks/Python.framework/Versions/7.2/lib/python2.7/site-packages/numpy/linalg/linalg.py", line 143, in _fastCopyAndTranspose
cast_arrays = cast_arrays + (_fastCT(a),)
MemoryError

[1]+  Exit 1                  python rand_diag.py

如何分配更多内存以及如何并行运行（我有 4 个内核）？

【问题讨论】：

这不是others 推荐的。
完全不清楚将矩阵乘以它的逆矩阵如何帮助“找到我的机器达到其极限的矩阵大小”。你真正想做什么？您是指内存限制、浮点舍入错误变得无法容忍的限制、执行时间限制还是其他什么？如果你想要一个单位矩阵，请使用numpy.identity(n)。
获取True 在教学上非常没有吸引力，因为它会导致学生将有限精度浮点运算与无限精度实数运算混淆。学习计算机逼近线性代数的学生需要学会不断意识到这种差异，并预测其后果。
我承认你的观点都非常有效。 @EricPostpischil：目前，我在运行 size=10000 矩阵的反转时观察到 100% 的 CPU 负载。内存仍然相当可用，如果我知道如何使用，我会将我机器的所有四个 CPU 分配给脚本，并为其提供所需的所有内存。有很多我想在这里更清楚地弄清楚......
@PatriciaShanahan：我目前最感兴趣，我们称之为规范线性代数，所以我只需要一些说明性的计算来了解正在发生的事情。计算细节稍后会引起关注，但不能无人看管。

标签： python numpy floating-point linear-algebra

【解决方案1】：

最后，你可以用

来四舍五入你的答案

m = np.round(m, decimals=10)

或检查它们是否有很大不同：

np.abs(A*A.I - i).mean() < 1e-10

如果你想杀死微小的数字。

我会用numpy.matrix 类来实现它。

import numpy

size = 100
A = numpy.matrix(numpy.random.randint(0,10,(size,)*2))
E = numpy.eye(size)

print A * A.I
print np.abs(A * A.I - E).mean() < 1e-10

【讨论】：

使用矩阵类的原因是什么？
它是有限的，但对于简单的二维线性代数来说更易于使用和阅读，例如，请参阅我的编辑。
或者，@TMOTTM，至少使用数组而不是循环遍历它。如果您愿意，您可以将我在示例中所做的所有操作都用于数组而不是矩阵（np.identity 给出了np.eye 的数组形式）当然除了最后的矩阵代数，但您可以在那里使用np.linalg，和你一样。

【解决方案2】：

同意已经提出的大部分观点。但是，我建议不要查看单个非对角线元素，而是取它们的 rms 总和；这在某种意义上反映了由于计算不完善而泄漏到非对角项中的“能量”。然后，如果您将此 RMS 数除以对角线项的总和，您将获得逆运算效果的指标。例如以下代码：

import numpy
import matplotlib.pyplot as plt
from numpy import mean, sqrt
N = 1000
R = numpy.zeros(N)

for size in range(50,N,50):

  A = numpy.zeros((size, size))
  E = numpy.zeros((size, size))

  for i in range(size):
      for j in range(size):
          A[i][j]+=numpy.random.randint(10)
          if i == j:
              E[i][j]=1

  A_inv = numpy.linalg.linalg.inv(A)
  D = numpy.dot(A, A_inv) - E
  S = sqrt(mean(D**2))
  R[size] = S/size
  print "size: ", size, "; rms is ",  S/size

plt.plot(range(50,N,50), R[range(50, N, 50)])
plt.ylabel('RMS fraction')
plt.show()

表明 rms 误差非常稳定，数组的大小一直到 950x950 的大小（它确实减慢了一点......）。但是，它从来都不是“精确的”，并且存在一些异常值（大概是当矩阵更接近奇异时 - 这可能发生在随机矩阵中。）

示例图（每次运行它都会看起来有点不同）：

【讨论】：

+1 鼓舞人心的实验，顺便说一句，你是量子化学家吗？
谢谢。实际上，我是一名实验物理学家/工程师。
我猜……“能量”一定来自某个地方……:)

【解决方案3】：

虽然获得True 具有教学吸引力，但它也会脱离浮点计算的现实。

在处理浮点数时，不仅要为不精确的结果做好准备，还要为出现的各种其他数值问题做好准备。

在您的特定情况下，为确保 A * inv(A) 足够接近单位矩阵，您可以计算 matrix norm 的 numpy.dot(A, A_inv) - E 并确保它足够小。

附带说明，您不必使用循环来填充A 和E。相反，您可以使用

A = numpy.random.randint(0, 10, (size,size))
E = numpy.eye(size)

【讨论】：

请注意，您建议的矩阵范数类似于我在答案中计算的 RMS 值 - 不同之处在于我然后除以矩阵的大小以获得相对比例的感觉。 +1 非常有用的参考资料！
+1 表示矩阵人口的捷径，我正在等待有人指出这一点。

【解决方案4】：

您的问题可以简化为常见的浮点比较问题。比较此类数组的正确方法是：

EPS = 1e-8  # for example
(np.abs(numpy.dot(A, A_inv) - E) < EPS).all()

【讨论】：