计算机是怎么求解线性方程的(矩阵乘和求逆)

上回我们说到，高斯老哥用消元法解线性方程，大致步骤呢就是给系数矩阵消元，运气好点呢直接整出上三角系数矩阵，得到方程组的唯一解，运气不行呢，消着消着发现整不出上三角，这时就得再讨论方程是有多解还是无解。 这里所说的"运气"呢其实可以根据行列式啊， A x = 0 Ax=0 Ax=0是否有解啊判断得到，具体操作可以看看我聊消元法的那一篇文章。 但是，高斯消元法存在一个问题，就是它是给人做的，比如给第一行乘个倍数加到另一行，或者将矩阵的两行交换个位置，这些我们手写计算当然没什么问题，但是你让计算机做它就不能忍了，计算机没那么多判断的能力，它就想要一种统一的计算模式和计算次数。ok，那我们就来聊一聊计算机是怎么做这些操作的，顺便我们能够知道为什么有很多科学家对矩阵乘法性能的提高有着孜孜不倦的追求。 我们先以一种抽象的方式来看看矩阵和列向量相乘，比如

[ 1 3 4 2 7 5 4 8 9 ] ∗ [ x y z ] = x ∗ [ 1 2 4 ] + y ∗ [ 3 7 8 ] + z ∗ [ 4 5 9 ] \begin{bmatrix} 1 & 3 & 4\\ 2 & 7 & 5 \\ 4 & 8 & 9 \end{bmatrix}* \begin{bmatrix} x \\ y \\ z \end{bmatrix}= x*\begin{bmatrix} 1 \\ 2 \\ 4 \end{bmatrix}+ y*\begin{bmatrix} 3 \\ 7 \\ 8 \end{bmatrix}+ z*\begin{bmatrix} 4 \\ 5 \\ 9 \end{bmatrix} ⎣⎡​124​378​459​⎦⎤​∗⎣⎡​xyz​⎦⎤​=x∗⎣⎡​124​⎦⎤​+y∗⎣⎡​378​⎦⎤​+z∗⎣⎡​459​⎦⎤​

根据第一节所说，系数方阵的三个列向量为三维空间中的三个独立向量，用它们可以组成三维空间的任一向量， [ x , y , z ] [x,y,z] [x,y,z]分别为这三个向量的扩展倍数，假如我想知道一倍的 [ 1 , 2 , 4 ] [1,2,4] [1,2,4]与三倍的 [ 3 , 7 , 8 ] [3,7,8] [3,7,8]的和向量是什么，只需要给 [ x , y , z ] [x,y,z] [x,y,z]赋值 [ 1 , 3 , 0 ] [1,3,0] [1,3,0]就可以了，即

[ 1 3 4 2 7 5 4 8 9 ]  ⁣ ∗  ⁣ [ 1 3 0 ]  ⁣ =  ⁣ 1  ⁣ ∗  ⁣ [ 1 2 4 ]  ⁣ +  ⁣ 3  ⁣ ∗  ⁣ [ 3 7 8 ]  ⁣ +  ⁣ 0  ⁣ ∗  ⁣ [ 4 5 9 ]  ⁣ =  ⁣ [ 10 23 28 ] \begin{bmatrix} 1 & 3 & 4\\ 2 & 7 & 5 \\ 4 & 8 & 9 \end{bmatrix}\!*\! \begin{bmatrix} 1 \\ 3 \\ 0 \end{bmatrix}\!=\! 1\!*\!\begin{bmatrix} 1 \\ 2 \\ 4 \end{bmatrix}\!+\! 3\!*\!\begin{bmatrix} 3 \\ 7 \\ 8 \end{bmatrix}\!+\! 0\!*\!\begin{bmatrix} 4 \\ 5 \\ 9 \end{bmatrix}\!=\!\begin{bmatrix} 10 \\ 23 \\ 28 \end{bmatrix} ⎣⎡​124​378​459​⎦⎤​∗⎣⎡​130​⎦⎤​=1∗⎣⎡​124​⎦⎤​+3∗⎣⎡​378​⎦⎤​+0∗⎣⎡​459​⎦⎤​=⎣⎡​102328​⎦⎤​

那么如果我想知道 [ 1 , 2 , 4 ] [1,2,4] [1,2,4]向量与 负三倍的 [ 1 , 2 , 4 ] [1,2,4] [1,2,4]与一倍的 [ 3 , 7 , 8 ] [3,7,8] [3,7,8]的和向量以及 负四倍的 [ 1 , 2 , 4 ] [1,2,4] [1,2,4]与一倍的 [ 4 , 5 , 9 ] [4,5,9] [4,5,9]的和向量 的组合是什么，就可以运算 [ 1 3 4 2 7 5 4 8 9 ]  ⁣ ∗  ⁣ [ 1 0 0 ] =    1  ⁣ ∗  ⁣ [ 1 2 4 ]  ⁣ +  ⁣ 0  ⁣ ∗  ⁣ [ 3 7 8 ]  ⁣ +  ⁣ 0  ⁣ ∗  ⁣ [ 4 5 9 ]  ⁣ ⁣ =  ⁣ [ 1 2 4 ] \begin{bmatrix} 1 & 3 & 4\\ 2 & 7 & 5 \\ 4 & 8 & 9 \end{bmatrix}\!*\! \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}\enspace= \enspace\;1\!*\!\begin{bmatrix} 1 \\ 2 \\ 4 \end{bmatrix}\!+\! 0\!*\!\begin{bmatrix} 3 \\ 7 \\ 8 \end{bmatrix}\!+\! 0\!*\!\begin{bmatrix} 4 \\ 5 \\ 9 \end{bmatrix}\!\!=\!\begin{bmatrix} 1 \\ 2\\ 4 \end{bmatrix} ⎣⎡​124​378​459​⎦⎤​∗⎣⎡​100​⎦⎤​=1∗⎣⎡​124​⎦⎤​+0∗⎣⎡​378​⎦⎤​+0∗⎣⎡​459​⎦⎤​=⎣⎡​124​⎦⎤​ [ 1 3 4 2 7 5 4 8 9 ]  ⁣ ∗  ⁣ [ − 3 1 0 ]  ⁣ =  ⁣ − 3  ⁣ ∗  ⁣ [ 1 2 4 ]  ⁣ +  ⁣ 1  ⁣ ∗  ⁣ [ 3 7 8 ]  ⁣ +  ⁣ 0  ⁣ ∗  ⁣ [ 4 5 9 ]  ⁣ ⁣ =  ⁣ [ 0 1 − 4 ] \begin{bmatrix} 1 & 3 & 4\\ 2 & 7 & 5 \\ 4 & 8 & 9 \end{bmatrix}\!*\! \begin{bmatrix} -3 \\ 1 \\ 0 \end{bmatrix}\!=\! -3\!*\!\begin{bmatrix} 1 \\ 2 \\ 4 \end{bmatrix}\!+\! 1\!*\!\begin{bmatrix} 3 \\ 7 \\ 8 \end{bmatrix}\!+\! 0\!*\!\begin{bmatrix} 4 \\ 5 \\ 9 \end{bmatrix}\!\!=\!\begin{bmatrix} 0 \\ 1\\ -4 \end{bmatrix} ⎣⎡​124​378​459​⎦⎤​∗⎣⎡​−310​⎦⎤​=−3∗⎣⎡​124​⎦⎤​+1∗⎣⎡​378​⎦⎤​+0∗⎣⎡​459​⎦⎤​=⎣⎡​01−4​⎦⎤​ [ 1 3 4 2 7 5 4 8 9 ]  ⁣ ∗  ⁣ [ − 4 0 1 ]  ⁣ =  ⁣ − 4  ⁣ ∗  ⁣ [ 1 2 4 ]  ⁣ +  ⁣ 0  ⁣ ∗  ⁣ [ 3 7 8 ]  ⁣ +  ⁣ 1  ⁣ ∗  ⁣ [ 4 5 9 ]  ⁣ ⁣ =  ⁣ [ 0 − 3 − 7 ] \begin{bmatrix} 1 & 3 & 4\\ 2 & 7 & 5 \\ 4 & 8 & 9 \end{bmatrix}\!*\! \begin{bmatrix} -4 \\ 0 \\ 1 \end{bmatrix}\!=\! -4\!*\!\begin{bmatrix} 1 \\ 2 \\ 4 \end{bmatrix}\!+\! 0\!*\!\begin{bmatrix} 3 \\ 7 \\ 8 \end{bmatrix}\!+\! 1\!*\!\begin{bmatrix} 4 \\ 5 \\ 9 \end{bmatrix}\!\!=\!\begin{bmatrix} 0 \\ -3\\ -7 \end{bmatrix} ⎣⎡​124​378​459​⎦⎤​∗⎣⎡​−401​⎦⎤​=−4∗⎣⎡​124​⎦⎤​+0∗⎣⎡​378​⎦⎤​+1∗⎣⎡​459​⎦⎤​=⎣⎡​0−3−7​⎦⎤​

将计算过程及结果整理一下就可以得到

[ 1 3 4 2 7 5 4 8 9 ]  ⁣ ∗  ⁣ [ 1 − 3 − 4 0 1 0 0 0 1 ]  ⁣ =  ⁣ [ 1 0 0 2 1 − 3 4 − 4 − 7 ] \begin{bmatrix} 1 & 3 & 4\\ 2 & 7 & 5 \\ 4 & 8 & 9 \end{bmatrix}\!*\! \begin{bmatrix} 1&-3&-4\\ 0&1&0\\ 0&0&1 \end{bmatrix}\!=\! \begin{bmatrix} 1 & 0 & 0\\ 2 & 1 & -3 \\ 4 & -4 & -7 \end{bmatrix} ⎣⎡​124​378​459​⎦⎤​∗⎣⎡​100​−310​−401​⎦⎤​=⎣⎡​124​01−4​0−3−7​⎦⎤​

此时我们会惊奇得发现居然消元了，第一行居然消元了！只不过在列方向，高斯消元法需要的是在行方向做消元，那如果我对行向量能做一些求和之类的操作问题不就解决了吗。而这你也不需要担心，因为有数学大佬已经帮你总结好了，即

[ x , y , z ] [ 1 3 4 2 7 5 4 8 9 ] = x [ 1 , 3 , 4 ] + y [ 2 , 7 , 5 ] + z [ 4 , 8 , 9 ] \begin{bmatrix} x ,y ,z \end{bmatrix} \begin{bmatrix} 1 & 3 & 4\\ 2 & 7 & 5 \\ 4 & 8 & 9 \end{bmatrix}=x[1,3,4]+y[2,7,5]+z[4,8,9] [x,y,z​]⎣⎡​124​378​459​⎦⎤​=x[1,3,4]+y[2,7,5]+z[4,8,9]

根据消元步骤先根据pivot1给23行消元，即

[ 1 0 0 − 2 1 0 − 4 0 1 ] [ 1 3 4 2 7 5 4 8 9 ] = [ 1 3 4 0 1 − 3 0 − 4 − 7 ] \begin{bmatrix} 1&0&0\\ -2&1&0\\ -4&0&1 \end{bmatrix} \begin{bmatrix} 1 & 3 & 4\\ 2 & 7 & 5 \\ 4 & 8 & 9 \end{bmatrix}= \begin{bmatrix} 1 & 3 & 4\\ 0 & 1 & -3 \\ 0 & -4 & -7 \end{bmatrix} ⎣⎡​1−2−4​010​001​⎦⎤​⎣⎡​124​378​459​⎦⎤​=⎣⎡​100​31−4​4−3−7​⎦⎤​

再根据pivot2给第三行消元，即

[ 1 0 0 0 1 0 0 4 1 ] [ 1 3 4 0 1 − 3 0 − 4 − 7 ] = [ 1 3 4 0 1 − 3 0 0 − 19 ] \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 4 & 1 \end{bmatrix} \begin{bmatrix} 1 & 3 & 4\\ 0 & 1 & -3 \\ 0 & -4 & -7 \end{bmatrix}= \begin{bmatrix} 1 & 3 & 4\\ 0 & 1 & -3 \\ 0 & 0 & -19 \end{bmatrix} ⎣⎡​100​014​001​⎦⎤​⎣⎡​100​31−4​4−3−7​⎦⎤​=⎣⎡​100​310​4−3−19​⎦⎤​

再将最后一行变为1

[ 1 0 0 0 1 0 0 0 − 1 19 ] [ 1 3 4 0 1 − 3 0 0 − 19 ] = [ 1 3 4 0 1 − 3 0 0 1 ] \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & -\dfrac{1}{19} \end{bmatrix} \begin{bmatrix} 1 & 3 & 4\\ 0 & 1 & -3 \\ 0 & 0 & -19 \end{bmatrix}= \begin{bmatrix} 1 & 3 & 4\\ 0 & 1 & -3 \\ 0 & 0 & 1 \end{bmatrix} ⎣⎢⎡​100​010​00−191​​⎦⎥⎤​⎣⎡​100​310​4−3−19​⎦⎤​=⎣⎡​100​310​4−31​⎦⎤​

以上我们将消元的过程变为了矩阵乘法，我猜这就是计算机做消元的方法，当然，在消元过程中有时要用到两行交换位置的情况，我们也可以利用矩阵乘得到，比如交换矩阵的前两行

[ 0 1 0 1 0 0 0 0 1 ] [ 1 3 4 2 7 5 4 8 9 ] = [ 2 7 5 1 3 4 4 8 9 ] \begin{bmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} 1 & 3 & 4\\ 2 & 7 & 5 \\ 4 & 8 & 9 \end{bmatrix}= \begin{bmatrix} 2 & 7 & 5 \\ 1 & 3 & 4\\ 4 & 8 & 9 \end{bmatrix} ⎣⎡​010​100​001​⎦⎤​⎣⎡​124​378​459​⎦⎤​=⎣⎡​214​738​549​⎦⎤​

我们将上述给原矩阵做变换的矩阵叫做初等变换矩阵，初等变换矩阵的逆矩阵非常好求，逆矩阵的定义为 A A − 1 = E AA^{-1}=E AA−1=E，而初等变换矩阵正是 E E E经过了简单的运算得到的。比如将第一行乘4并加到第二行的变换矩阵，只要从第二行中减去第一行的4倍就变回了单位阵，也就得到了逆矩阵。，即

A = [ 1 0 0 4 1 0 0 0 1 ] ⇒ A − 1 = [ 1 0 0 − 4 1 0 0 0 1 ] A=\begin{bmatrix} 1 & 0 & 0 \\ 4 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}⇒ A^{-1}=\begin{bmatrix} 1 & 0 & 0 \\ -4 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} A=⎣⎡​140​010​001​⎦⎤​⇒A−1=⎣⎡​1−40​010​001​⎦⎤​

而交换两行的变换阵的逆矩阵是它本身

那么更一般的矩阵求逆怎么做呢(已知矩阵可逆)，要不怎么说数学家牛逼呢，Gauss-Jordan法使求逆完全可以通过消元来解决，就是这二位 他们表示，假如你想求一个逆矩阵啊，先给原矩阵旁边整一个单位阵，然后给第一个矩阵开始消元，消元过程中单位阵做同样的变换，当原来的矩阵变为单位阵时，单位阵就是原矩阵的逆矩阵，首先是Gauss的消元法创造上三角矩阵

[ 1 3   ∣   1 0 2 7   ∣   0 1 ] ⇒ [ 1 3   ∣   1 0 0 1   ∣   − 2 1 ] \begin{bmatrix} 1 & 3 \,|\,1&0\\ 2 & 7 \,|\,0&1 \end{bmatrix}⇒ \begin{bmatrix} 1 & 3 \,|\,1&0\\ 0 & 1 \,|\,-2&1 \end{bmatrix} [12​3∣17∣0​01​]⇒[10​3∣11∣−2​01​]

然后jordan说你把这个上三角再倒着消元成单位阵就得到逆矩阵啦

[ 1 3   ∣   1 0 0 1   ∣   − 2 1 ] ⇒ [ 1 0   ∣   7 − 3 0 1   ∣   − 2 1 ] \begin{bmatrix} 1 & 3 \,|\,1&0\\ 0 & 1 \,|\,-2&1 \end{bmatrix}⇒ \begin{bmatrix} 1 & 0 \,|\,7&-3\\ 0 & 1 \,|\,-2&1 \end{bmatrix} [10​3∣11∣−2​01​]⇒[10​0∣71∣−2​−31​]

不放心的话拿到matlab里验证一下就ok了。 我们可以看到，化简求逆都是经过矩阵乘来运算，这种较为统一的运算过程可是把计算机爽死了，我们说的解方程啊，矩阵求逆啊到这里都变成了矩阵的乘法运算，而矩阵乘不仅仅只有这一点应用，在空间中一个向量的放大，缩小，平移，旋转都是通过乘矩阵的方式来实现，所以在矩阵乘这个计算上一点点的优化都可能使软件性能发生巨大变化，这也就回答了文章开始的问题，为什么有那么多人不断追求矩阵乘运算性能的提高。微小的改变却能导致巨大的性能变化，这可能就是算法之美吧。