用Python求矩阵的广义逆

对于两个方阵$A,B$，若$AB=E$，且$E$为单位阵，则$A,B$互逆，可记作A=B−1,B=A−1A=B^{-1}, B=A^{-1}A=B−1,B=A−1。

inv

在numpy和scipy中，均提供了求逆函数，分别是numpy.linalg.inv和scipy.lingalg.inv，下面举个例子看一下

import numpy as np
import numpy.linalg as nl
import scipy.linalg as slA = np.random.rand(3,3)

二者求逆的结果如下，是完全相同的

>>> nl.inv(A)
array([[  0.56940181,   7.07580181,  -2.80811204],[ -0.4241864 ,  -0.45826635,   1.51146259],[  1.80687147, -13.30035928,   2.91842045]])
>>> sl.inv(A)
array([[  0.56940181,   7.07580181,  -2.80811204],[ -0.4241864 ,  -0.45826635,   1.51146259],[  1.80687147, -13.30035928,   2.91842045]])

下面测试以下二者的计算时间

from timeit import timeit
A = np.random.rand(500,500)
timeit(lambda : nl.inv(A), number=100)
# 2.689510200000001
timeit(lambda : sl.inv(A), number=100)
# 0.5326764000000139

可见scipy.linalg.inv要比numpy更快。

穆尔-彭罗斯广义逆

有了逆这个概念，那么对于求解$Ax=b$这样的方程组就变得简单了，$x$可写为x=A−1bx=A^{-1}bx=A−1b。

但并不是所有矩阵都有逆的，为了让x=A−1bx=A^{-1}bx=A−1b这种形式拥有更广的适用范围，就必须拓展矩阵的逆的概念，此即广义逆。

对于复矩阵$A$而言，若存在复矩阵$G$，满足

1. $AGA=A$
2. $GAG=G$
3. (AG)H=AG(AG)^H=AG(AG)H=AG
4. (GA)H=GA(GA)^H=GA(GA)H=GA

其中∗H*^H∗H表示共轭转置，则称$G$为$A$的穆尔-彭罗斯广义逆，也叫加号逆，记作G=A+G=A^+G=A+。

scipy.linalg提供了两个求加号逆的矩阵，分别是pinv和pinvh，分别用于实矩阵和复矩阵。

pinv

在inv中，除了待求矩阵a之外，还有两个参数：overwrite_a默认为False，若为True，则在计算过程中覆盖a；check_finite默认False，若为True则进行有限性检查。

在pinv, pinvh中，除了a和check_finite之外，还有下列参数：

• atol, rtol 用于控制精度
• return_rank 若为True，则返回矩阵的有效秩

下面以pinv为例，简单介绍一下函数的使用方法。

A = np.random.rand(3,4)
# sl.inv(A)会报错
G = sl.pinv(A)

接下来看一下$A$和$AGA$是否相等

>>> A
array([[0.06135389, 0.53545499, 0.52057907, 0.48222946],[0.65828574, 0.42574002, 0.09521576, 0.62040737],[0.86266928, 0.51081037, 0.11678121, 0.74462639]])
>>> A@G@A
array([[0.06135389, 0.53545499, 0.52057907, 0.48222946],[0.65828574, 0.42574002, 0.09521576, 0.62040737],[0.86266928, 0.51081037, 0.11678121, 0.74462639]])

可见二者是相等的。