「算法笔记」矩阵求逆

发表于 | | 阅读次数:
字数统计: 655 字
矩阵求逆和整数逆元比较类似，是通过高斯消元的思路实现的。

思路

求 $n$ 阶矩阵 $A$ 的逆矩阵，也就是 $A^{-1}$，$\frac{1}{A}$，求出一个 $n$ 阶矩阵 $B$ 使得 $AB=E$（其中 $E​$ 表示单位矩阵）。

先回顾一下矩阵的三种初等行变换（具体内容详见 「算法笔记」矩阵入门）：

  1. 交换某两行。
  2. 将某一行的所有元素乘上 $k$（其中 $k\neq 0$）。
  3. 将某一行的所有元素乘上 $k$ 后加到另一行上去。

假设我们通过 $P_1,P_2,P_3,\dots,P_k$ 等初等矩阵分别左乘 $A$，即 $P_1P_2P_3\cdots P_kA=PA$，其中的 $P$ 为 $P_1$ 到 $P_k$ 的乘积。我们想让 $A$ 变为 $P$ 非常简单，只需要使用高斯消元先变成上三角矩阵，然后变成对角矩阵。

考虑如何求出这个 $P$，我们首先有 $PA=E$，$PE=P$，如果我们同时维护两个矩阵 $A$ 和 $B$ 并使得 $B$ 一开始等于 $E$，在对 $A$ 进行初等行变换变为 $E$ 的过程中也对 $B$ 进行相同的初等行变换。由于初等行变换等价于被对应的初等矩阵左乘，所以当 $A$ 变成 $P$ 后，$B$ 也就从 $E$ 变成了 $P$。

对于无解情况，如果在高斯消元的过程中发现 $A$ 无法变成 $E$ 就说明无解。

时间复杂度：$O(n^3)$

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
 
#include <cstdio>
#include <cmath>
#include <algorithm>

const int N=305;
const int mod=1e9+7;
int n,m,a[N][N<<1];

int pow(int x,int p) {
    int ret=1;
    for(;p;p>>=1,x=1LL*x*x%mod) if(p&1) ret=1LL*ret*x%mod;
    return ret;
}
bool Gauss(int n,int m) {
    for(int i=1;i<=n;++i) {
        int p=i;
        for(int k=i+1;k<=n;++k) if(std::abs(a[k][i])>std::abs(a[p][i])) p=k;
        if(i!=p) std::swap(a[i],a[p]);
        if(!a[i][i]) return 0;
        int d=pow(a[i][i],mod-2);
        for(int j=i;j<=m;++j) a[i][j]=1LL*a[i][j]*d%mod;
        for(int k=1;k<=n;++k) {
            if(i==k) continue;
            int d=a[k][i];
            for(int j=i;j<=m;++j) a[k][j]=(a[k][j]-1LL*a[i][j]*d%mod+mod)%mod;
        }
    }
    return 1;
}
int main() {
    scanf("%d",&n),m=n+n;
    for(int i=1;i<=n;++i) {
        for(int j=1;j<=n;++j) scanf("%d",&a[i][j]);
        a[i][n+i]=1;
    }
    if(!Gauss(n,m)) return puts("No Solution"),0;
    for(int i=1;i<=n;++i) {
        for(int j=n+1;j<=m;++j) printf("%d%c",a[i][j]," \n"[j==m]);
    }
    return 0;
}

习题

0%