一、伪逆概述

        并非所有矩阵都有逆矩阵。逆用于求解方程组，在某些情况下，方程组没有解，因此逆不存在。然而，找到一个最小化误差的近似值是有意义的。例如，使用伪逆找到一组数据点的最佳拟合线。

        Moore-Penrose 伪逆是SVD的直接应用。

        正如我们在之前篇章中看到的，矩阵A的逆矩阵可用于求解方程Ax=b：

        

         但是在方程组有0个或多个解的情况下，无法找到逆，也无法求解方程，这时候需要伪逆上场了。

        伪逆记作​​​​​​​，并，并最小化

         以下公式可用于求伪逆：

        

        （1）U、D 和 V 分别是 A 的左奇异向量、奇异值和右奇异向量。

        （2）是A的伪逆，是D的伪逆。

        （3）我们知道D是对角矩阵，因此D+可以通过取D的非零值的倒数来计算。

 二、例1：伪逆计算

         让我们看看如何计算它。 我们将创建一个非方阵 A，计算其奇异值分解及其伪逆。

        

A = np.array([[7, 2], [3, 4], [5, 3]])
U, D, V = np.linalg.svd(A)

D_plus = np.zeros((A.shape[0], A.shape[1])).T
D_plus[:D.shape[0], :D.shape[0]] = np.linalg.inv(np.diag(D))

A_plus = V.T.dot(D_plus).dot(U.T)

        求得A_plus

array([[ 0.16666667, -0.10606061,  0.03030303],
       [-0.16666667,  0.28787879,  0.06060606]])

          我们现在可以使用 Numpy 的 pinv() 函数检查伪逆是否正确：

np.linalg.pinv(A)

        求得

array([[ 0.16666667, -0.10606061,  0.03030303],
       [-0.16666667,  0.28787879,  0.06060606]])

        数据一致，我们现在可以检查它是否真的是 A 的近似倒数。

        我们已知并且

A_plus.dot(A)

         求得下面的矩阵，可以看出近似单位矩阵，对吧。

array([[  1.00000000e+00,   2.63677968e-16],
       [  5.55111512e-17,   1.00000000e+00]])

        就是说，也就是

        计算伪逆的另一种方法是使用以下公式：

        

A_plus_1 = np.linalg.inv(A.T.dot(A)).dot(A.T)

        对于不同矩阵结果不如SVD方法精准。Numpy的函数pinv()就是使用SVD。

三、例2：使用伪逆求解超定线性方程组

        对于这个例子，我们将考虑这组三个具有两个未知数的方程：

        

        让我们可视化

x1 = np.linspace(-5, 5, 1000)
x2_1 = -2*x1 + 2
x2_2 = 4*x1 + 8
x2_3 = -1*x1 + 2

plt.plot(x1, x2_1)
plt.plot(x1, x2_2)
plt.plot(x1, x2_3)
plt.xlim(-2., 1)
plt.ylim(1, 5)
plt.show()

         可视化效果如下，三条线没有交与一点，就是没有解。伪逆从最小二乘误差的角度求解系统：它找到最小化误差的解。

         我们用矩阵形式表示

        

        即

        我们现在将计算 A 的伪逆：

A = np.array([[-2, -1], [4, -1], [-1, -1]])
A_plus = np.linalg.pinv(A)
A_plus

         求得

array([[-0.11290323,  0.17741935, -0.06451613],
       [-0.37096774, -0.27419355, -0.35483871]])

         取四位小数即，

         我们可以用它来找到 x：、

b = np.array([[-2], [-8], [-2]])
res = A_plus.dot(b)
res

         求得

array([[-1.06451613],
       [ 3.64516129]])

         即x坐标

        让我们绘制到图上看一下

plt.plot(x1, x2_1)
plt.plot(x1, x2_2)
plt.plot(x1, x2_3)
plt.xlim(-2., 1)
plt.ylim(1, 5)

plt.scatter(res[0], res[1])

plt.show()

         可视化效果如下

         也许您会期望该点位于三角形的重心处。但情况并非如此，因为方程的缩放方式不同。

四、例3：直线拟合

        此方法还可用于将线拟合到一组点。设有以下数据点：

         我们现有这组 x 和 y，我们要寻找最小化误差的线 y=mx+b。 可以将误差评估为拟合和实际数据点之间的差异之和。

        我们可以用矩阵方程表示数据点：

        

         请注意，这里的矩阵 A 表示系数的值。 1 列对应于截距（没有它，拟合将有跨越原点的约束）。 它给出了以下方程组：

        

        我们有一组方程 mx+b=y。 这些用于返回截距参数。 例如，在对应于第一个点的第一个方程中，我们有很好的 x=0 和 y=2。 这可能会令人迷惑，因为这里的向量 x 对应于系数。 这是因为我们正在寻找直线的系数而不是 x 和 y 未知数。

         所以我们将构造这些矩阵并尝试使用伪逆来找到最小化误差（线与实际数据点之间的差异）的线的方程。

        让我们创建 A 和 b 的矩阵并计算A的伪逆：

A = np.array([[0, 1], [1, 1], [2, 1], [3, 1], [3, 1], [4, 1]])
b = np.array([[2], [4], [0], [2], [5], [3]])
A_plus = np.linalg.pinv(A)

        求得伪逆如下

array([[ -2.00000000e-01,  -1.07692308e-01,  -1.53846154e-02,
          7.69230769e-02,   7.69230769e-02,   1.69230769e-01],
       [  6.00000000e-01,   4.00000000e-01,   2.00000000e-01,
          4.16333634e-17,   4.16333634e-17,  -2.00000000e-01]])

        带入公式x=A+b

coefs = A_plus.dot(b)

        求得拟合的参数。 斜率为 m=0.21538462，截距为 b=2.2。

array([[ 0.21538462],
       [ 2.2       ]])

        我们将绘制数据点和回归线：

x = np.linspace(-1, 5, 1000)
y = coefs[0]*x + coefs[1]

plt.plot(A[:, 0], b, '*')
plt.plot(x, y)
plt.xlim(-1., 6)
plt.ylim(-0.5, 5.5)

plt.show()

        使用numpy库进行核验，使用polyfit函数对A和b进行拟合，得到如下结果，说明计算无误。

np.polyfit(A[:, 0], b, 1)

array([[ 0.21538462],
       [ 2.2       ]])

         还有很多计算库或者在线拟合网页可以验证上面的计算，可以自行尝试。

五、例4：直线拟合

        要查看具有更多数据点的过程，我们可以生成数据。

        我们将生成一个包含 100 个具有随机 x 值和伪随机 y 值的点的列向量（参见下面的 reshape()）。 Numpy.random 包中的函数 seed() 用于冻结随机化并能够重现结果：

np.random.seed(123)
x = 5*np.random.rand(100)
y = 2*x + 1 + np.random.randn(100)

x = x.reshape(100, 1)
y = y.reshape(100, 1)

        我们将通过添加一列 1 从 x 创建矩阵 A，就像我们在示例 3 中一样。

A = np.hstack((x, np.ones(np.shape(x))))
A[:10]

        得到

array([[ 3.48234593,  1.        ],
       [ 1.43069667,  1.        ],
       [ 1.13425727,  1.        ],
       [ 2.75657385,  1.        ],
       [ 3.59734485,  1.        ],
       [ 2.1155323 ,  1.        ],
       [ 4.90382099,  1.        ],
       [ 3.42414869,  1.        ],
       [ 2.40465951,  1.        ],
       [ 1.96058759,  1.        ]])

        我们现在可以找到 A 的伪逆并计算回归线的系数：

A_plus = np.linalg.pinv(A)
coefs = A_plus.dot(y)
coefs

        求得

array([[ 1.9461907 ],
       [ 1.16994745]])

        我们将点和直线可视化

x_line = np.linspace(0, 5, 1000)
y_line = coefs[0]*x_line + coefs[1]

plt.plot(x, y, '*')
plt.plot(x_line, y_line)
plt.show()

六、其它参考

        您可以看到伪逆对这类问题非常有用！另外关于多项式拟合或拟合平面可以参考以下文章。

Opencv学习笔记 - 多项式求解和拟合_bashendixie5的博客-CSDN博客_opencv 多项式拟合曲线的函数类型有:双曲线型、对数型、指数型、多项式型等。对于最常见的曲线, 其特征接近于多项式型, 所以这里的曲线拟合问题就变成了多项式回归求解.一般情况下,对实际的点数据进行求解是无解的, 所以我们需要引入最小二乘法来求解。https://blog.csdn.net/bashendixie5/article/details/115176277Opencv学习笔记 - 使用最小二乘法拟合平面_bashendixie5的博客-CSDN博客_opencv 平面拟合https://blog.csdn.net/u012719076/article/details/113124887https://blog.csdn.net/liumangmao1314/article/details/54179526本文主要验证了博客上的最小二乘法拟合平面的。与 用matlab拟合出来的平面计算的点到直线的距离是一样的，而且系数也是一样的。说明了本方法的可行性。matlab中公式为z = c + ax +byoepncv中公式为Ax+By+Cz=D 将openc..https://blog.csdn.net/bashendixie5/article/details/115413982