3.5 矩阵广义逆

定义

  矩阵的一号逆，是对矩阵逆的推广。对于线性方程组$A\mathbf{x}=\mathbf{b}$来说，如果$A$可逆，那么它的解就是$x=A^{-1}b$。然后对于不可逆的方程，可不可以找到一个矩阵$G$，使得$G\mathbf{b}$是方程的解呢？对于相容方程组，答案是肯定的，这个矩阵$G$就是矩阵的一号逆$1-inverse$，也叫减号逆，或广义逆generalized inverse，记号为$A^{-}$。它的定义是：
$$A{A}^{-}A=A$$
  符合这个定义的所有矩阵都是矩阵的减号逆，对于可逆矩阵来说，减号逆只有一个，对于不可逆矩阵来说，减号逆存在多个。

算法

  根据定义是很难求出减号逆的。求减号逆的办法是将矩阵先扩充变成一个方阵，变成下面这个样子：
$$\left(\begin{array}{cc}{A}_{m\times n}& E_{m\times m}\\ E_{n\times n}& 0\end{array}\right)$$
  通过填充把矩阵变成了$(m+n)\times (m+n)$的矩阵，再进行初等行变换与初等列变换，将左上角变成单位矩阵$E_r$,其中$r$是矩阵的秩，这样这个扩充的矩阵就变成这个样子：
$$\left(\begin{array}{cc}{E}_r& P\\ Q& 0\end{array}\right)$$
  最后的减号逆就是一个矩阵乘法就出来了：
$$A^{-}=Q{\left(\begin{array}{cc}{E}_r& \ast \\ \ast & \ast \end{array}\right)}_{n\times m}P$$
  这里的$\ast$表示放什么数字都可以了。为了计算方便，可以全部取0。

举例

  以下是矩阵：
$$A=\left(\begin{array}{ccc}1& 1& 2\\ 2& 2& 1\\ 3& 3& 3\end{array}\right)$$
  扩展成大矩阵：
$$\left(\begin{array}{cccccc}1& 1& 2& 1& 0& 0\\ 2& 2& 1& 0& 1& 0\\ 3& 3& 3& 0& 0& 1\\ 1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\end{array}\right)$$
  先初等行变换：
$$\left(\begin{array}{cccccc}1& 1& 2& 1& 0& 0\\ 0& 0& 1& 0.667& -0.333& 0\\ 0& 0& 0& -1& -1& 1\\ 1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\end{array}\right)$$
  再进行初等列变换：
$$\left(\begin{array}{cccccc}1& 0& 0& 1& 0& 0\\ 0& 1& 0& 0.667& -0.333& 0\\ 0& 0& 0& -1& -1& 1\\ 1& -2& -1& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\end{array}\right)$$
  所以得到了$P,Q$,以下就是：
$$\begin{gathered} P=\left(\begin{array}{ccc}1& 0& 0\\ 0.667& -0.333& 0\\ -1& -1& 1\end{array}\right) \\ Q=\left(\begin{array}{ccc}1& -2& -1\\ 0& 0& 1\\ 0& 1& 0\end{array}\right) \end{gathered}$$
  $E_r$剩余部位全为0的减号逆为：
$$A^{-}=Q{\left(\begin{array}{cc}{E}_r& 0\\ 0& 0\end{array}\right)}_{n\times m}P=\left(\begin{array}{ccc}-0.333& 0.667& 0\\ 0& 0& 0\\ 0.667& -0.333& 0\end{array}\right)$$

齐次方程组通解

  通过矩阵的减号逆，可以获取齐次方程组的所有解。虽然之前说过求齐次方程组通解的办法。但是这里给出另外一种办法，这种办法就是利用矩阵的减号逆，齐次方程组$A\mathbf{x}=0$的通解是：
$$(E_n-A^{-}A)\mathbf{y},\mathbf{y}\in R^n$$
  公式中的$y$是任意变量，不过这种形式表示通解不太方便。因为还需要继续展开计算，才可以变成自由变量与基的线性组合。

非齐次方程组通解

  同样非齐次方程组的解也可以用减号逆表示，我们知道非齐次方程组的解就是特解加上齐次方程组的通解，所以对于$A\mathbf{x}=\mathbf{b}$，所有的解就是：
$$A^{-}\mathbf{b}+(E_n-A^{-}A)\mathbf{y},\mathbf{y}\in R^n$$
  同样，需要展开表示。

Python实现

  减号逆的实现，我用python实现了一下,只贴了部分代码:

   # 减号逆
    def one_inverse(self):

        form, echelons = self.to_p_q_form()
        # 再获取Q矩阵，P矩阵
        m = len(self.__vectors[0])
        n = len(self.__vectors)
        p_array = [[0 for _ in range(m)] for _ in range(m)]
        for column in range(n, m + n):
            for line in range(0, m):
                p_array[column - n][line] = form.__vectors[column][line]
        q_array = [[0 for _ in range(n)] for _ in range(n)]
        for column in range(0, n):
            for line in range(m, m + n):
                q_array[column][line - m] = form.__vectors[column][line]
        er_array = [[0 for _ in range(n)] for _ in range(m)]

        for i, echelon in enumerate(echelons):
            if echelon[0] >= n:
                break
            er_array[i][i] = 1

        q = Matrix(q_array)
        er = Matrix(er_array)
        p = Matrix(p_array)
        return (q * er )* p

    # pq型
    def to_p_q_form(self):
        m = len(self.__vectors[0])
        n = len(self.__vectors)

        new_matrix = self.expand_q_r_form(m, n)
        new_matrix = new_matrix.to_upper_triangle(last_line=m)

        # 变成阶梯形后，每一行除于自己的主元系数
        echelons = new_matrix.get_echelons()
        vectors = new_matrix.__vectors
        for line, echelon in enumerate(echelons):
            column = echelon[0]
            if column >= n:
                break
            matrix_utils.row_times(vectors, line, 1 / vectors[column][line])

        # 循环阶梯进行初等列变换,将右面都变成0
        for line, echelon in enumerate(echelons):
            # 把右边全变成0
            main_column = echelon[0]
            pivot = vectors[main_column][line]
            for column in range(main_column + 1, n):
                line_pivot = -vectors[column][line] / pivot
                matrix_utils.column_times_add(vectors, main_column, line_pivot, column)


        # 再进行初等列变换,将主元列移动到前面
        for line, echelon in enumerate(echelons):

            main_column = echelon[0]
            if main_column >= n:
                break
            matrix_utils.column_swap(vectors, main_column, line)
        return new_matrix, echelons

    def expand_q_r_form(self, m, n):
        # 组建一个新矩阵
        vectors = [[0 for _ in range(m + n)] for _ in range(m + n)]
        # 先将矩阵复制进去
        for i, vector in enumerate(self.__vectors):
            for j, e in enumerate(vector):
                vectors[i][j] = e
        # 右边放一个
        for column in range(n, m + n):
            vectors[column][column - n] = 1
        # 下面放一个
        for row in range(m, m + n):
            vectors[row - m][row] = 1
        # 行变换到阶梯形
        new_matrix = Matrix(vectors)
        return new_matrix