如何从纽比.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)

运行代码会产生

^{pr2}$

很明显，结果是一个单位矩阵，但并不精确，因此print numpy.dot(A, A_inv) == E显然给出了{}。我这样做是为了练习线性代数，并试图找出矩阵的大小，我的机器达到了它的极限。得到一个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

同意已经提出的大部分观点。然而，我建议不要看个别的非对角元素，而是取它们的均方根和；这在某种意义上反映了由于计算不完善而泄漏到非对角项中的“能量”。如果你把这个均方根数除以对角线项的和，你就得到了一个度量逆运算效果的度量。例如，以下代码：

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的情况下非常稳定（它确实减慢了一点…）。然而，它从来都不是“精确的”，并且存在一些异常值（可能是当矩阵更接近奇异时-随机矩阵也可能发生这种情况）

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

Answer 2

虽然得到True在教学上很有吸引力，但它也将脱离浮点计算的实际情况。在

在处理浮点时，我们不仅要为不精确的结果做好准备，而且还要为出现的所有其他数值问题做好准备。在

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

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

作为补充说明，您不必使用循环来填充A和{}。相反，你可以用

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

Answer 3

最后，你可以用

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

或者看看他们是否有很大的不同：

^{pr2}$

如果你想消灭这些微小的数字。在

---

我将使用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