从矩阵乘法到高斯消元

· · 个人记录

前置知识：一点点的数学基础。

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]} 配套练习题： [矩阵乘法](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项

最基本的方法就是O(n)递推，显然这题会超时。现在我就要介绍这一章节的题目，矩阵快速幂了。

说到矩阵快速幂，首先要学习这一题，如果没有掌握的话，先去学习一下。

我们的题目要求其实就是：

\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初始化成这样就行。

至于乘法的重定义，就是我们上一章节讲的【矩阵乘法】的内容，这里便不再赘述。

矩阵快速幂的模板：

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 【模板】矩阵快速幂

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 【模板】乘法逆元

（首先理解常数的逆元，更有助于理解矩阵求逆。）

先定义矩阵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}

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}

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}

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}

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}矩阵的逆矩阵公式：

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}是指第一个方程，第一个未知数的系数，以此类推。）

将未知数的系数和方程的值写成矩阵：

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 【模板】高斯消元法

对于这一题而言，我们可以选择高斯-约旦消元法。

约旦消元法大致思路如下：

  1. 选择一个尚未被选过的未知数，选择一个包含这个主元的方程。
  2. 将这个方程主元的系数化为1。
  3. 通过加减消元，消掉其它方程中的这个未知数。
  4. 重复以上步骤，直到把每一行都变成只有一项有系数。

相比普通的高斯消元，约旦消元法的代码更加简单，且没有回带的过程。

这一题的代码：

#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]线性方程组

对于这一题而言，比模板多了一个判断无解的情况，那么用约旦消元法不太好判断无解，便可以用朴素的高斯消元法。

高斯消元法的主要步骤与高斯-约旦消元法的步骤差不多，就是少了一个把与这个未知数有关的所有的系数都交换位置的操作。再消完元后，高斯-约旦消元法可以直接求解，而高斯消元法需要回带。

正因为这一点不同，所以高斯-约旦消元法消完后是只有主对角线上有非0数的形状，可以省去回带过程，类似于：

a_{11} & 0 & …… & 0 & b_1\\ 0 & a_{22} & …… & 0 & b_2\\ …… &……&……&……\\ 0& 0 & …… & a_{nn} &b_n \end{bmatrix}

所以说可以直接将非0的那一项等于右边的b数组中对应的项。

但正因如此，高斯-约旦消元法不太好判断无解。

而高斯消元法消完之后是一个类似于倒三角形的形状。

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，这个方程组无解。

（注意，这一题要先判断无解，再判断多解。）

对于考场上的题目来说，个人建议如果能用的话，尽量用高斯-约旦消元法，因为他的码量相比还是比较小的。

这一题的代码：

#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]外星千足虫

对于这一题而言，我们只需要把消元的时候的除法换成异或即可。且因为是异或，可以用bitset来优化，缩小时间复杂度。

代码（仅供参考）：

#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] 1&2 \end{bmatrix}

那么根据这个原理，我们可以得出以下求出矩阵的逆的步骤：

  1. 先把矩阵A和I放在一起。
  2. 将A变成单位矩阵。
  3. 消完元后的I即是A^{-1}。

（其实求矩阵的逆有很多种方法，但是我tcl觉得这种是最好懂的，就写了这种。）

代码：

#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。