Linear Algebra(5):初等矩阵和矩阵的逆

线性代数

Elementary Matrix

线性系统角度求解矩阵的逆

首先回归一下什么是矩阵的逆：矩阵中$AB=BA=I$ ，则称B是A的逆矩阵，记做$B=A^{-1}$ .从这个定义可以发现只有方阵才有逆矩阵.那么现在假设方阵A存在逆矩阵，那么如何求解这个逆矩阵呢？

比如矩阵A为下面的形式：

$$A=\begin{pmatrix} 1 & 2 \\3& 4\end{pmatrix} \quad A.A^{-1}=I=\begin{pmatrix} 1 & 0 \\0& 1\end{pmatrix} \quad$$

此时求解$A^{-1}$ 一种比较直观的想法就是先假设出来：

$$A^{-1} = \begin{pmatrix}x_{11} &x_{12}\\ x_{21}&x_{22} \end{pmatrix} \\$$

接着依据矩阵乘法的定义可得：

$$\begin{pmatrix}1 &2\\ 3&4 \end{pmatrix}.\begin{pmatrix}x_{11} &x_{12}\\ x_{21}&x_{22} \end{pmatrix} = \begin{pmatrix}1 &0\\ 0&1 \end{pmatrix} \\$$

最终可得到有四个未知数的四个方程组：

$$\begin{cases} x_{11}+2x_{21}=1\\ 3x_{11}+4x_{21}=0\\ x_{12}+2x_{22}=0\\ 3x_{12}+4x_{22}=1 \end{cases}$$

显然求解矩阵的逆本质可以规约成一个求解线性系统的问题， 四个未知数的四个方程，这个解要么有唯一解或无数解或无解.当有唯一解这个结果就是矩阵相应的逆，如果无解就代表矩阵不可逆，无数解呢？其实对于这样的一个线性系统是不可能有无数解的.

如下图所示，把这个矩阵乘法拆成两部分，左边对应$(1 \quad0)^T$ ,右边对应$(0\quad1)^T$, 问题就变成了分别求解左右两个方程组.左边只含有$x_{11}、x_{21}$ 两个未知数，而右边只含有$x_{12}、x_{22}$ 两个未知数，使用Gauss-Jordan消元法将两边的系数矩阵化为行最简形，其中的*号就是对应的解.

同时可以发现这个过程中，左右两个方程组的系数矩阵是相同的，区别只在于右边的结果向量.左边对应$(1 \quad0)^T$ ,右边对应$(0\quad1)^T$.而Gauss-Jordan消元法就是对系数矩阵进行操作，因此可以将两个方程组合一块(如上图)，直接对上面的增广矩阵化简，如果系数矩阵能化为行最简形，那么右侧就是相应的解$x_{11}、x_{12}、x_{21}、x_{22}$ .

接下来就具体求解一下：

$$\begin{pmatrix} 1 & 2 &1&0\\3& 4&0&1\end{pmatrix} \quad \rightarrow \begin{pmatrix} 1 & 2 &1&0\\0& -2&-3&1\end{pmatrix} \quad \rightarrow \begin{pmatrix} 1 & 2 &1&0\\0& 1&1.5&-0.5\end{pmatrix} \quad \rightarrow \begin{pmatrix} 1 & 0 &-2&1\\0& 1&1.5&-0.5\end{pmatrix} \quad \\$$

因此:

$$A^{-1}=\begin{pmatrix}-2&1\\1.5&-0.5\end{pmatrix} \quad \\$$

观察下面这个增广矩阵显然是不可能有无数解的，因为不存在全0行. 但是可能是无解的，此时矩阵A没有逆矩阵.

$$\left( \begin{array}{cc|cc} a & b & 1 & 0 \\ c & d & 0 & 1 \end{array} \right)$$

也与前面讲到矩阵逆的性质对应即：对于矩阵A，如果存在逆矩阵B，则B唯一. 那么何时无解？显然是当系数矩阵化为行最简形式时存在0行.

代码实现：

求解矩阵A的逆的过程中增广矩阵右侧不再是一个列向量b，而可能是和A同型的一个单位矩阵,因此这里对构造函数进行了一定的改造：

from .Matrix import Matrix
from .Vector import Vector
from ._global import is_zero

class LinearSystem:
    def __init__(self, A, b):
        assert A.row_num() == len(b),\
        "row number of A must be equal to  the length of b"
        self._m = A.row_num()
        self._n = A.col_num()

        if isinstance(b, Vector):
            self.Ab = [Vector(A.row_vector(i).underlying_list() + [b[i]]) for i in range(self._m)]
        if isinstance(b, Matrix):
            self.Ab = [Vector(A.row_vector(i).underlying_list() + b.row_vector(i).underlying_list()) for i in range(self._m)]  
        self.pivots = []
    
    def _max_row(self, index_i, index_j, n):
        best, ret= self.Ab[index_i][index_j], index_i
        for i in range(index_i+1, n):
            if self.Ab[i][index_j] > best:
                best, ret = self.Ab[i][index_j], i
        return ret

    def _forward(self):
        i, k = 0, 0
        while i < self._m and k < self._n: 
            # 看Ab[i][k]是否为主元
            max_row = self._max_row(i, k, self._m)
            self.Ab[i], self.Ab[max_row] = self.Ab[max_row], self.Ab[i]

            if is_zero(self.Ab[i][k]) :
                k += 1
            else:
                self.Ab[i] = self.Ab[i]/self.Ab[i][k]
                for j in range(i+1, self._m):
                    self.Ab[j] = self.Ab[j] - self.Ab[j][k]*self.Ab[i]
                self.pivots.append(k)
                i+=1
    
    def _backward(self):
        n = len(self.pivots)
        for i in range(n-1, -1, -1): 
            k = self.pivots[i]
            for j in range(i-1, -1, -1):
                self.Ab[j] = self.Ab[j] - self.Ab[j][k]*self.Ab[i]

    def gauss_jordan_elimination(self):

        self._forward()
        self._backward()

        for i in range(len(self.pivots), self._m):
            if not is_zero(self.Ab[i][-1]):
                return False
        return True

    def fancy_print(self):
        for i in range(self._m):
            print(" ".join(str(self.Ab[i][j])for j in range(self._n)),end=" ")
            print("|", self.Ab[i][-1])

def inv(A):

    if A.row_num()!=A.col_num():
        return None

    n = A.row_num()
    ls = LinearSystem(A, Matrix.identity(n))
    ls.gauss_jordan_elimination()

    if not ls.gauss_jordan_elimination:
        return None
    
    invA = [[row[i] for i in range(n, 2*n)] for row in ls.Ab]

    return Matrix(invA)

测试代码：

from LA.LinearSystem import LinearSystem, inv
from LA.Matrix import Matrix
from LA.Vector import Vector


if __name__ == "__main__":
    A = Matrix([[1,2],[3,4]])
    invA = inv(A)
    print(invA)

可以发现这个结果就是上面我们求解的结果.

初等矩阵

我们解决一个线性系统的过程其实就是对表示整个线性系统的矩阵进行下面的三种矩阵操作：

• 矩阵的某一行乘以一个常数
• 矩阵的一行加(减)另一行
• 交换矩阵的两行

这三种基本的矩阵操作还是在矩阵中针对特定的元素去操作。但同时矩阵本身可以表示变换，像上篇中在二维空间中图形变换的例子就是，只需用一个变换矩阵T去乘以另外一个可以表示点集的矩阵，将这些点转换成了另外的一些点.这样的操作是在矩阵的级别上的操作，我们要想进行这种变换，只需要对表示点集的矩阵在左乘一个变换矩阵就好了，而不是在点集矩阵中逐一的元素操作.

那么上面求解线性系统的三种操作是否也可以用矩阵表示？这个问题可以从单位矩阵出发.

首先单位矩阵满足下面的式子：

$$\begin{pmatrix} 1 &0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\\ \end{pmatrix} . \begin{pmatrix} a &b & c\\ d & e & f\\ g & h & i\\ \end{pmatrix}= \begin{pmatrix} a &b & c\\ d & e & f\\ g & h & i\\ \end{pmatrix}$$

对单位矩阵稍加改造就能得到我们想要的初等矩阵.

1.矩阵的某一行乘以一个常数

$$\begin{pmatrix} k &0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\\ \end{pmatrix} . \begin{pmatrix} a &b & c\\ d & e & f\\ g & h & i\\ \end{pmatrix}= \begin{pmatrix} ka &kb & kc\\ d & e & f\\ g & h & i\\ \end{pmatrix}$$

2.矩阵的一行加(减)另一行的若干倍

$$\begin{pmatrix} 1 &0 & 0\\ 0 & 1 & 0\\ 1 & 0 & 1\\ \end{pmatrix} . \begin{pmatrix} a &b & c\\ d & e & f\\ g & h & i\\ \end{pmatrix}= \begin{pmatrix} a &b & c\\ d & e & f\\ a+g & b+h & c+i\\ \end{pmatrix}$$

$$\begin{pmatrix} 1 &0 & -1\\ 0 & 1 & 0\\ 0 & 0 & 1\\ \end{pmatrix} . \begin{pmatrix} a &b & c\\ d & e & f\\ g & h & i\\ \end{pmatrix}= \begin{pmatrix} a-g &b-h & c-i\\ d & e & f\\ g & h & i\\ \end{pmatrix}$$

$$\begin{pmatrix} 1 &0 & 0\\ 0 & 1 & -p\\ 0 & 0 & 1\\ \end{pmatrix} . \begin{pmatrix} a &b & c\\ d & e & f\\ g & h & i\\ \end{pmatrix}= \begin{pmatrix} a &b & c\\ d-pg & e-ph & f-pi\\ g & h & i\\ \end{pmatrix}$$

3.交换矩阵的两行

$$\begin{pmatrix} 1 &0 & 0\\ 0 & 0 & 1\\ 0 & 1 & 0\\ \end{pmatrix} . \begin{pmatrix} a &b & c\\ d & e & f\\ g & h & i\\ \end{pmatrix}= \begin{pmatrix} a &b & c\\ g & h & i\\ d & e & f\\ \end{pmatrix}$$

总结三种矩阵，如果我们要让结果矩阵进行其中变换的一种，只需要对单位矩阵进行对应的变换就好了. 这三中基本操作有一个专门的名词即矩阵的初等变换.

矩阵的初等变换都可以表示成上面其中的一种变换矩阵，这些对应矩阵的就叫做初等矩阵. 初等矩阵即对单位矩阵进行一次初等变换得到的结果矩阵.通常记作$E$ .

而前面的Gauss-Jordan消元法把矩阵化为行最简形式的过程就是寻找一系列初等矩阵E，使得：

$$E_p...E_3.E_2.E_1.A=rref(A)$$

初等矩阵角度的可逆性

因为初等变换是可逆的，所以初等矩阵是可逆的。

以三类初等变换为例：

• 矩阵的某一行乘以一个常数

$$E_1=\begin{pmatrix} k &0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\\ \end{pmatrix} \quad E_2= \begin{pmatrix} \frac{1}{k} &0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\\ \end{pmatrix}$$

如果想要逆变这个操作，只需要在同一行乘以一个1/k.显然E1和E2互为逆.

• 矩阵的一行加(减)另一行的若干倍

$$E_1=\begin{pmatrix} 1 &0 & 0\\ 0 & 1 & -p\\ 0 & 0 & 1\\ \end{pmatrix} \quad E_2= \begin{pmatrix} 1 &0 & 0\\ 0 & 1 & p\\ 0 & 0 & 1\\ \end{pmatrix}$$

E1为第二行减去第三行的p倍，显然逆操作就是E2将第二行加上第三行的p倍.

• 交换矩阵两行

$$E_1=\begin{pmatrix} 1 &0 & 0\\ 0 & 0 & 1\\ 0 & 1 & 0\\ \end{pmatrix} \quad E_2= \begin{pmatrix} 1 &0 & 0\\ 0 & 0 & 1\\ 0 & 1 & 0\\ \end{pmatrix}$$

E1为交换第二行和第三行，显然逆操作就是E2依旧交换第二行和第三行. 可以发现E1和E2相同.

所以所有的初等矩阵都是可逆的.

那么一个一般的矩阵如果它可逆相应的如何得到逆矩阵呢？

如果一个矩阵A可逆，存在一系列E,使得：

$$E_p....E_3.E_2.E_1.A=I\\ E_p....E_3.E_2.E_1.A.A^{-1}=I.A^{-1}\\ E_p....E_3.E_2.E_1.I=A^{-1}\\ \left( \begin{array}{c|c} A&I \\ \end{array} \right) \rightarrow \left( \begin{array}{c|c} I&A^{-1} \\ \end{array} \right)$$

矩阵的逆为什么重要

首先矩阵的逆有一个非常重要的计算功能，比如下面的线性系统：

$$Ax=b\\ x=A^{-1}.b$$

其次矩阵的逆和很多重要的命题连接在了一起。即对于方阵A，下面这些这些命题都是等价的：

$$矩阵A可逆 \quad A是非奇异矩阵\iff 线性系统Ax=0 只有唯一解x=0 \iff rref(A)=I\iff A可以表示成一系列初等矩阵的乘积 \\ \iff Ax=b只有唯一解$$

矩阵的LU分解

我们在小学时期学过数字的分解，任意一个整数分解成若干个素数相乘：

$$30=2*3*5$$

同理一个矩阵也可以分解为几个矩阵乘积的形式，不同的矩阵分解方式有不同的目的，这里矩阵的LU分解是以提高计算效率为目的.

矩阵的LU分解就是将一个矩阵A分解为：

$$A=L.U$$

其中L代表Lower Triangle Matrix下三角矩阵，而U代表Upper Triangle Matrix上三角矩阵.

什么是上三角矩阵(左)和下三角矩阵（中）?以方阵为例：

$$\begin{pmatrix} * &0 & 0 &0\\ * &* & 0 &0\\ * &* & * &0\\ * &* & * &*\\ \end{pmatrix} \quad \begin{pmatrix} * &* & * &*\\ 0 &* & * &*\\ 0 &0 & * &*\\ 0 &0 & 0 &*\\ \end{pmatrix} \quad \begin{pmatrix} 1 &0 & 0 &0\\ * &1 & 0 &0\\ * &* & 1 &0\\ * &* & * &1\\ \end{pmatrix} \quad$$

矩阵的LU分解就是将一个矩阵分成这样两个矩阵的乘积，更具体来说其中的L是单位下三角矩阵(右)，主对角线值都为1.

回忆高斯消元的过程：

$$\begin{pmatrix} 1 &2 & 4\\ 3 & 7 & 2\\ 2 & 3 & 3\\ \end{pmatrix} \quad \rightarrow \begin{pmatrix} 1 &2 & 4\\ 0 & 1 & -10\\ 0 & -1 & -5\\ \end{pmatrix} \quad \rightarrow \begin{pmatrix} 1 &2 & 4\\ 0 & 1 & -10\\ 0 & 0 & -15\\ \end{pmatrix} \quad$$

高斯消元法的过程就是通过一系列初等变换把一个矩阵变成了上三角矩阵.

$$E_p...E_3.E_2.E_1.A=U\\ A=E_1^{-1}.E_2^{-1}.E_3^{-1}...E_p^{-1}.U$$

一切顺利的话此时的L就是一个下三角矩阵就有：

$$L=E_1^{-1}.E_2^{-1}.E_3^{-1}...E_p^{-1}$$

因此并不是所有的矩阵都可以进行LU分解，而是有条件的.

以下面的方阵为例:

• 高斯消元过程

$$\begin{pmatrix} 1 &2 & 3\\ 4 & 5 & 6\\ 3 & -3 & 5\\ \end{pmatrix}\quad \rightarrow \begin{pmatrix} 1 &2 & 3\\ 0 & -3 & -6\\ 3 & -3 & 5\\ \end{pmatrix}\quad \rightarrow \begin{pmatrix} 1 &2 & 3\\ 0 & -3 & -6\\ 0 & -9 & -4\\ \end{pmatrix}\quad \rightarrow \begin{pmatrix} 1 &2 & 3\\ 0 & -3 & -6\\ 0 & 0 & 14\\ \end{pmatrix}\quad$$

• 对应的初等矩阵逆矩阵过程

$$\begin{pmatrix} 1 &0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\\ \end{pmatrix}\quad \rightarrow \begin{pmatrix} 1 &0 & 0\\ 4 & 1 & 0\\ 0 & 0 & 1\\ \end{pmatrix} \quad \rightarrow \begin{pmatrix} 1 &0 & 0\\ 4 & 1 & 0\\ 3 & 0 & 1\\ \end{pmatrix} \quad \rightarrow \begin{pmatrix} 1 &0 & 0\\ 4 & 1 & 0\\ 3 & 3 & 1\\ \end{pmatrix} \quad$$

可以发现高斯消元法的过程是把下对角线化为0，而对应初等变换逆过来就会得到一个下三角矩阵.

因此:

$$A=\begin{pmatrix} 1 &2 & 3\\ 4 & 5 & 6\\ 3 & -3 & 5\\ \end{pmatrix}\quad L =\begin{pmatrix} 1 &0 & 0\\ 4 & 1 & 0\\ 3 & 3 & 1\\ \end{pmatrix}\quad U= \begin{pmatrix} 1 &2 & 3\\ 0 & -3 & -6\\ 0 & 0 & 14\\ \end{pmatrix}\quad$$

从上面过程可以发现矩阵可以进行LU分解的条件：对A进行高斯消元的过程中不需要交换两行的位置.

为什么要进行矩阵的LU分解？

在解$Ax=b$ 的过程中，对A进行LU分解需要复杂度$O(0.5n^3)$ ，此时变为$L.U.x=b$ 在设定$U.x=y$ ，那么$L.y=b$ 只需要$O(0.5n^2)$ 求出y. 然后$Ux=y$ 也是需要$O(0.5n^2)$ 求出x. 对比求解矩阵的逆需要$O(2n^3)$ .

LU分解代码实现

代码实现：

from .Matrix import Matrix
from .Vector import Vector
from ._global import is_zero

def lu(matrix):

    assert matrix.row_num() == matrix.col_num(),\
    "matrix must be a square matrix"

    n = matrix.row_num()
    A = [matrix.row_vector(i) for i in range(n)]
    L = [[1.0 if i==j else 0.0 for i in range(n)]for j in range(n)]

    for i in range(n):
        if is_zero(A[i][i]):
            return None, None
        else:
            for j in range(i+1, n):
                p = A[j][i]/A[i][i]
                A[j] = A[j] - p*A[i]
                L[j][i] = p
    return Matrix(L), Matrix([A[i].underlying_list() for i in range(n)])

测试代码,测试案例就是上面的例子：

from LA.LU import lu
from LA.Matrix import Matrix

if __name__=="__main__":
    m1 = Matrix([[1,2,3],[4,5,6],[3,-3,5]])
    L, U = lu(m1)
    print(L)
    print(U)

可以发现打印结果与我们上面的结果一致，说明代码正确.

LU分解还可以用于非方阵：

$$A=\begin{pmatrix} 1 &0 & 0 &0\\ * &1 & 0 &0\\ * &* & 1 &0\\ * &* & * &1\\ \end{pmatrix} \quad \begin{pmatrix} * &* & * &* &* & * \\ 0 &* & * &* &* & * \\ 0 &0 & * &* &* & * \\ 0 &0 & 0 &* &* & * \\ \end{pmatrix} \quad$$

此时的分解过程和方阵的分解过程是完全一样的.

但如果A是一个高大于宽的矩阵会稍微不一样：

$$A=\begin{pmatrix} 1 &0 & 0 &0\\ * &1 & 0 &0\\ * &* & 1 &0\\ * &* & * &1\\ * &* & * &* \\ * &* & * &* \\ \end{pmatrix} \quad \begin{pmatrix} * &* & * &* \\ 0 &* & * &* \\ 0 &0 & * &* \\ 0 &0 & 0 &* \\ \end{pmatrix} \quad$$