从矩阵乘法到高斯消元

#### 前置知识：一点点的数学基础。 ## 0. 矩阵的定义 一个$m\times{n}$的矩阵就是$m\times{n}$个数排成$m$行$n$列的一个数阵。 ~~没错这个就是用来水段落的，~~ 下面开始正文。 ---- ## 1. 矩阵乘法 一个矩阵的行数必须等于一个矩阵的列数，这两个矩阵相乘才有意义。 矩阵$A$规模是$n×m$，矩阵$B$规模是$m×p$，现在需要你求$A\times{B}$。 公式为： $$C[i][j]=\sum_{k=1}^{m} A[i][k] \times{B[k][j]}$$ $C$矩阵就是$A\times{B}$的结果。 配套练习题： [矩阵乘法](https://www.luogu.com.cn/problem/U152117) ---- ## 2. 矩阵快速幂 首先我们要先了解矩阵乘法的性质： 矩阵乘法满足结合律。$(A\times{B})\times{C}=A\times({B}\times{C})$ 矩阵乘法**不满足**交换律。$A\times{B}\neq B\times{A}$ （其实仔细想想，也就很好理解。） 然后我们放一道例题： [Fibonacci第n项](https://www.luogu.com.cn/problem/U152118) 最基本的方法就是$O(n)$递推，显然这题会超时。现在我就要介绍这一章节的题目，矩阵快速幂了。 说到矩阵快速幂，首先要学习[这一题](https://www.luogu.com.cn/problem/P1226)，如果没有掌握的话，先去学习一下。 我们的题目要求其实就是： $$\begin{bmatrix}\ Fn,F(n+1)\end{bmatrix}->\begin{bmatrix}\ F(n+1),Fn+F(n+1)\end{bmatrix}$$ 那么我们可不可以让$\begin{bmatrix}\ Fn,F(n+1)\end{bmatrix}$乘一个矩阵得到$\begin{bmatrix}\ F(n+1),Fn+F(n+1)\end{bmatrix}$呢？ 显然是可以的，容易得出这个矩阵是：$\begin{bmatrix} 1&1 \\ 1&0 \end{bmatrix}$ 设$\begin{bmatrix}\ F1,F2\end{bmatrix}$为$A$，$\begin{bmatrix} 1&1 \\ 1&0 \end{bmatrix}$为$B$ 那么我们就可以将问题转换为$A\times{B^n}$ 这里出现了$B^n$，就可以用矩阵快速幂。我们可以推算，矩阵快速幂的时间复杂度其实是$O(k^3\times{log n})$的（$k$为矩阵的边长），符合这一题的数据范围。所以矩阵快速幂在这种$n$特别大，而矩阵的边长又特别小的类型的题目中非常适用。 讲了这么久怎么用矩阵快速幂，现在来说一说怎么用代码实现矩阵快速幂。 像平常的快速幂一样，我们要用一个$t$来得出答案。但是我们该如何初始化这个$t$呢？在常数快速幂的时候，它是$1$，但是在矩阵中，他显然不能是$1$这个数，而应该是一个矩阵。这时候便要引进一种主对角线上都是$1$，其他元素都是$0$的矩阵，这种矩阵乘其他矩阵的乘积都是原来那个矩阵的。所以我们就将$t$初始化成这样就行。 至于乘法的重定义，就是我们上一章节讲的【矩阵乘法】的内容，这里便不再赘述。 矩阵快速幂的模板： ```cpp struct node{ ll a[3][3]; } node quickmulti(node a,node b){ node tmp; memset(tmp.a,0,sizeof(tmp.a)); for(int i=1;i<=2;i++) for(int j=1;j<=2;j++) for(int k=1;k<=2;k++) tmp.a[i][j]=(tmp.a[i][j]+a.a[i][k]*b.a[j][k]%m)%m; return tmp; } node quickpow(ll n){ node tmp; memset(tmp.a,0,sizeof(tmp.a)); for(int i=1;i<=2;i++) tmp.a[i][i]=1; while(n){ if(n&1) tmp=quickmulti(f,tmp); f=quickmulti(f,f); n>>=1; } return tmp; } ``` （这里的$2$就是矩阵的边长，在不同的题目中会有不同的要求。） 配套练习题： [P3390 【模板】矩阵快速幂](https://www.luogu.com.cn/problem/P3390) ---- ## 3. 矩阵求逆 现在这里就要用到我们学过的**真·一点点**数学基础了。 小学的时候，我们会遇到类似于这样的方程： $$ax=b$$ 我们通常的方法是将左右都乘上$\frac{1}{a}$，然后方程就会变为： $$ax\times{\frac{1}{a}}=b\times{\frac{1}{a}}$$ $$x=\frac{b}{a}$$ 这样做的原理是$a\times{\frac{1}{a}}=1$，左边只剩下来一个$x$，我们就可以将$x$求出来啦。 但是如果对应到矩阵中，我们又该如何去做呢？ 还是原来那个例子， $$Ax=B$$ 我们该如何求出$\frac{1}{A}$呢？ （其实$\frac{1}{A}$就是$A^{-1}$，后文都以$A^{-1}$代替$\frac{1}{A}$） 这时候就要用到我们的主角，矩阵求逆了！ 前置知识：[P3811 【模板】乘法逆元](https://www.luogu.com.cn/problem/P3811) （首先理解常数的逆元，更有助于理解矩阵求逆。） 先定义矩阵$A$的逆矩阵为$A^{-1}$。 有矩阵方程： $$Ax=B$$ 左右乘$A^{-1}$ $$A\times{A^{-1}}\times{x}=B\times{A^{-1}}$$ $$x=B\times{A^{-1}}$$ 那么我们就可以解出$x$，但是怎么求出$A^{-1}$呢？我们就可以再次列出一个方程： $$A\times{A^{-1}}=I$$ （这里的$I$是前文说过的“矩阵中的$1$”，也就是单位矩阵，主对角线上都是$1$，其他元素都是$0$） ~~下文开始手玩矩阵求逆（，不感兴趣的同学可以直接跳到代码处（。~~ 写这个的原因主要是想让各位更好的理解原理，而不是只是单纯的会写代码。 ---- 如果矩阵$A$=$\begin{bmatrix} 2&1 \\ 5&3 \end{bmatrix}$，那么$I=\begin{bmatrix} 1&0 \\ 0&1 \end{bmatrix}$ 既然$A\times{A^{-1}}=I$，那么$A^{-1}$肯定也是一个$2\times{2}$的矩阵。 取$E1=\begin{bmatrix} \frac{1}{2}&0 \\ 0&1 \end{bmatrix}$ $E1\times{A}=\begin{bmatrix} 1&\frac{1}{2} \\ 5&3 \end{bmatrix}$ 这里的$E1$除了$[1,1]$上是$A[1][1]$的倒数，其他和单位矩阵格式相同。 获得矩阵的$[1,1]$正好和$I[1][1]$相等。 取$E2=\begin{bmatrix} 1&0 \\ -5&1 \end{bmatrix}$ $E2\times({E1\times{A}})=\begin{bmatrix} 1&\frac{1}{2} \\ 0&\frac{1}{2} \end{bmatrix}$ 这里的$E2$除了$[2][1]$上是$E1\times{A}$的$[2][1]$的相反数，其他和单位矩阵格式相同。 获得矩阵的$[2][1]$正好和$I[2][1]$相等。 取$E3=\begin{bmatrix} 1&0 \\ 0&2 \end{bmatrix}$ $E3\times[{E2\times({E1\times{A}})}]=\begin{bmatrix} 1&\frac{1}{2} \\ 0&1 \end{bmatrix}$ 这里的$E3$除了$[2][2]$上是$E2\times({E1\times{A}})$的$[2][2]$的倒数，其他与单位矩阵格式相同。 获得矩阵的$[2][2]$正好和$I[2][2]$相等。 取$E4=\begin{bmatrix} 1&-\frac{1}{2} \\ 0&1 \end{bmatrix}$ $E4\times{E3\times[{E2\times({E1\times{A}})}]}=\begin{bmatrix} 1&0 \\ 0&1 \end{bmatrix}$ 这里的$E4$除了$[1][2]$是$E3\times[{E2\times({E1\times{A}})}]$的$[1][2]$的相反数，其他和单位矩阵格式相同。 获得矩阵的$[1][2]$正好与$I[1][2]$相同。 那么$E4\times{E3\times[{E2\times({E1\times{A}})}]}=\begin{bmatrix} 1&0 \\ 0&1 \end{bmatrix}$，将$E4\times{E3\times{E2\times{E1}}}$，得到的就是$A^{-1}$ 验算一下，$A^{-1}=E4\times{E3\times{E2\times{E1}}}=\begin{bmatrix} 3&-1 \\ -5&2 \end{bmatrix}$ 发现$A^{-1}\times{A}$确实等于$I$，构造完成。 这里只是手玩了$2\times{2}$矩阵的逆矩阵，下面给出一个一般$2\times{2}$矩阵的逆矩阵公式： $\begin{bmatrix} a&b \\ c&d \end{bmatrix}$的逆矩阵$=\begin{bmatrix} \frac{d}{ad-bc}&\frac{-c}{ad-bc} \\ \frac{-b}{ad-bc}&\frac{a}{ad-bc} \end{bmatrix}$ （有兴趣的同学可以自行证明，这里就不再赘述）。 可以发现求逆矩阵的步骤比较复杂，但是我们可以通过计算机很简单的求出来。这只是对于$2\times{2}$的简单情况，更普遍的情况请看后文如何用代码实现。 --- ## 4. 线性方程求解和高斯消元 要想用代码实现矩阵求逆，就还要学习“线性方程求解与高斯消元”。 考虑以下线性方程： $ \left\{\begin{matrix} a_{11}x_1+&a_{12}x_2+……+&a_{1n}x_n =b_1\\ a_{21}x_1+&a_{22}x_2+……+&a_{2n}x_n =b_2\\ ……&……&……\\ a_{n1}x_1+&a_{n2}x_2+……+&a_{nn}x_n =b_n \end{matrix}\right.$ （$a_{11}$是指第一个方程，第一个未知数的系数，以此类推。） 将未知数的系数和方程的值写成矩阵： $\begin{bmatrix} a_{11}&a_{12} ……&a_{1n}&b_1 \\ a_{21}&a_{22} ……&a_{2n}&b_2 \\ ……&……&……\\ a_{n1}&a_{n2} ……&a_{nn}&b_n \end{bmatrix}$ 我们就可以用加减消元和代入消元来解出这个方程（大概在初一学一元二次方程组的时候就可以掌握了）。 例题：[P3389 【模板】高斯消元法](https://www.luogu.com.cn/problem/P3389) 对于这一题而言，我们可以选择高斯-约旦消元法。 约旦消元法大致思路如下： 1. 选择一个尚未被选过的未知数，选择一个包含这个主元的方程。 2. 将这个方程主元的系数化为1。 3. 通过加减消元，消掉其它方程中的这个未知数。 4. 重复以上步骤，直到把每一行都变成只有一项有系数。 相比普通的高斯消元，约旦消元法的代码更加简单，且没有回带的过程。 这一题的代码： ```cpp #include<bits/stdc++.h> using namespace std; #define eps 1e-9 double a[110][110]; int n; void Gauss(){ for(int i=1;i<=n;i++){ int mx=i; for(int j=i+1;j<=n;j++) if(fabs(a[j][i])>fabs(a[mx][i])) mx=j; for(int j=1;j<=n+1;j++) swap(a[i][j],a[mx][j]); if(fabs(a[i][i])<eps){ cout<<"No Solution"<<endl; return; } for(int j=1;j<=n;j++) if(j!=i){ double tmp=a[j][i]/a[i][i]; for(int k=i+1;k<=n+1;k++) a[j][k]-=a[i][k]*tmp; } } for(int i=1;i<=n;i++) printf("%.2lf\n",a[i][n+1]/a[i][i]); return; } int main(){ cin>>n; for(int i=1;i<=n;i++) for(int j=1;j<=n+1;j++) cin>>a[i][j]; Gauss(); return 0; } ``` 刚才介绍了高斯-约旦消元法，现在再来说一说朴素的高斯消元法。 例题：[P2455 [SDOI2006]线性方程组](https://www.luogu.com.cn/problem/P2455) 对于这一题而言，比模板多了一个判断无解的情况，那么用约旦消元法不太好判断无解，便可以用朴素的高斯消元法。 高斯消元法的主要步骤与高斯-约旦消元法的步骤差不多，就是少了一个把与这个未知数有关的所有的系数都交换位置的操作。再消完元后，高斯-约旦消元法可以直接求解，而高斯消元法需要回带。 正因为这一点不同，所以高斯-约旦消元法消完后是只有主对角线上有非0数的形状，可以省去回带过程，类似于： $\begin{bmatrix} a_{11} & 0 & …… & 0 & b_1\\ 0 & a_{22} & …… & 0 & b_2\\ …… &……&……&……\\ 0& 0 & …… & a_{nn} &b_n \end{bmatrix}$ 所以说可以直接将非$0$的那一项等于右边的$b$数组中对应的项。 但正因如此，高斯-约旦消元法不太好判断无解。 而高斯消元法消完之后是一个类似于倒三角形的形状。 $\begin{bmatrix} a_{11} & a_{12} & …… & a_{1n} & b_1\\ 0 & a_{22} & …… & a_{2n} & b_2\\ …… &……&……&……\\ 0& 0 & …… & a_{nn} &b_n \end{bmatrix}$ 显然最后一个方程可以直接求解，然后可以把最后一个方程的解代到倒数第二个方程中，这样倒数第二个方程也可以求解……依次类推，这就是回带的过程。 所以，高斯消元可以判断无解或多解，当这一个方程所有未知数的系数都为$0$，且对应的$b$也为$0$时，方程有多解。只要有一个方程的所有未知数系数都为$0$，且对应的$b$不为$0$，这个方程组无解。 （注意，这一题要先判断无解，再判断多解。） 对于考场上的题目来说，个人建议如果能用的话，尽量用高斯-约旦消元法，因为他的码量相比还是比较小的。 这一题的代码： ```cpp #include<bits/stdc++.h> using namespace std; const double eps=1e-9; double a[101][101],ans[101]; int n; void Gauss(){ for(int i=1;i<=n;i++){ int tab=i; for(int j=i+1;j<=n;j++) if(fabs(a[j][i])>fabs(a[tab][i])) tab=j; if(tab!=i) swap(a[tab],a[i]); if(fabs(a[i][i])>=eps) for(int j=1;j<=n;j++){ if(i==j) continue; double temp=a[j][i]/a[i][i]; for(int k=1;k<=n+1;k++) a[j][k]-=temp*a[i][k]; } } int s1=0,s2=0; for(int i=1;i<=n;i++){ int num=0; for(int j=1;j<=n+1;j++) if(a[i][j]==0) num++; else break; if(num==n+1) s1=1; if(num==n&&a[i][n+1]!=0) s2=1; } if(s2==1) {cout<<-1<<endl;return;} if(s1==1) {cout<<0<<endl;return;} for(int i=n;i>=1;i--){ ans[i]=a[i][n+1]/a[i][i]; for(int j=i;j>=1;j--) a[j][n+1]-=a[j][i]*ans[i]; } for(int i=1;i<=n;i++) printf("x%d=%.2lf\n",i,ans[i]); } int main(){ cin>>n; for(int i=1;i<=n;i++) for(int j=1;j<=n+1;j++) cin>>a[i][j]; Gauss(); return 0; } ``` 高斯消元求解异或方程组例题：[P2447 [SDOI2010]外星千足虫](https://www.luogu.com.cn/problem/P2447) 对于这一题而言，我们只需要把消元的时候的除法换成异或即可。且因为是异或，可以用$bitset$来优化，缩小时间复杂度。 代码（仅供参考）： ```cpp #include<bits/stdc++.h> using namespace std; bitset<2005>b[2005]; string s; int ans,n,m; void Gauss(){ for(int i=1;i<=n;i++){ int tab=i; while(tab<=m&&!b[tab][i]) ++tab; ans=max(ans,tab); if(ans>m) return; if(tab!=i) swap(b[i],b[tab]); for(int j=1;j<=m;j++) if(i!=j&&b[j][i]) b[j]^=b[i]; } } int main(){ cin>>n>>m; for(int i=1;i<=m;i++){ int x; cin>>s>>x; for(int j=1;j<=n;j++) b[i][j]=s[j-1]-'0'; b[i][n+1]=x; } Gauss(); if(ans>m){ cout<<"Cannot Determine"<<endl; return 0; } cout<<ans<<endl; for(int i=1;i<=n;i++) if(b[i][n+1]==0) cout<<"Earth"<<endl; else cout<<"?y7M#"<<endl; return 0; } ``` 相信如果你认真地看到了这里，一定对高斯消元有了一定的认识，现在我就来填之前挖的坑——矩阵求逆： 设矩阵$A$的逆矩阵为$A^{-1}$，单位矩阵为$I$ 那么$A^{-1}\times{[AI]}=[IA^{-1}]$ 这里的$[]$是指把两个矩阵放到一起，比如说 $A=[1],B=[2]$ $[AB]=\begin{bmatrix} 1&2 \end{bmatrix}$ 那么根据这个原理，我们可以得出以下求出矩阵的逆的步骤： 1. 先把矩阵$A$和$I$放在一起。 2. 将$A$变成单位矩阵。 3. 消完元后的$I$即是$A^{-1}$。 （其实求矩阵的逆有很多种方法，但是我~~tcl~~觉得这种是最好懂的，就写了这种。） 代码： ```cpp #include<bits/stdc++.h> using namespace std; #define ll long long #define mod 1000000007 ll a[501][1001],n; ll quickpow(ll n,ll m){ ll ans=1; while(m){ if(m&1) ans=ans*n%mod; n=n*n%mod; m>>=1; } return ans%mod; } void Gauss(){ for(int i=1;i<=n;i++){ int mx=i; for(int j=i+1;j<=n;j++) if(a[j][i]>a[mx][i]) mx=j; if(mx!=i) swap(a[i],a[mx]); if(!a[i][i]){ cout<<"No Solution"<<endl; return; } ll ny=quickpow(a[i][i],mod-2); for(int k=1;k<=n;k++){ if(k==i) continue; ll p=a[k][i]*ny%mod; for(int j=i;j<=n*2;j++) a[k][j]=((a[k][j]-p*a[i][j])%mod+mod)%mod; } for(int j=1;j<=n*2;j++) a[i][j]=(a[i][j]*ny%mod); } for(int i=1;i<=n;i++){ for(int j=n+1;j<=n*2;j++) cout<<a[i][j]<<" "; cout<<endl; } return; } int main(){ cin>>n; for(int i=1;i<=n;i++) for(int j=1;j<=n;j++){ cin>>a[i][j]; a[i][i+n]=1; } Gauss(); return 0; } ``` ## 5. 小结 矩阵在现实生活中有很多的作用，比如说乘法原理和加法原理（就是四年级奥数上的那个）都与矩阵有着密不可分的关系。 在数学上，矩阵远远不止我这篇文章的所介绍的作用。比如说：矩阵在概率上的应用，矩阵的几何性质等等等等…… 如果感兴趣的话可以自行BFS，如果本人水准达到这个程度的时候可能还会再写一篇关于矩阵的文章~~wtcl~~。 写文章不易，不喜勿喷。如果文章中有错误请私信告诉我，只要我看到，立刻就会更正。文章都是原创，没有转载任何文献。如有雷同纯属巧合，并请私信告诉我，我来修改一下。 求评论qwq求点赞qwq。