python - 如何从 numpy.dot(A, A_inv) 得到单位矩阵

Question

我准备了一个随机数矩阵，计算它的逆矩阵并将它与原始矩阵相乘。这在理论上给出了单位矩阵。我怎样才能让 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 个内核)？

Answer 1

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

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

我强烈推荐阅读 What Every Computer Scientist Should Know About Floating-Point Arithmetic .

在您的特定情况下，为确保 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)