矩阵求逆 —— 初等变换法（高斯-约旦消元）

一只萌新

·

2019-08-21 17:26:36

·

题解

\color{teal}{Blog}

P4783 【模板】矩阵求逆

题目描述

求一个N×N的矩阵的逆矩阵。答案对10^9+7取模。

1.逆矩阵的定义

假设 A 是一个方阵，如果存在一个矩阵 A^{-1}，使得

A^{-1}A=I

并且

AA^{-1}=I

那么，矩阵 A 就是可逆的，A^{-1} 称为 A 的逆矩阵

2.逆矩阵求法 —— 初等变换法（高斯-约旦消元）

0.高斯-约旦消元

详见

P3389 【模板】高斯消元法

题解部分

高斯约旦消元与高斯消元区别：

高斯消元 -> 消成上三角矩阵

高斯-约旦消元 -> 消成对角矩阵

约旦消元法的精度更好,代码更简单,没有回带的过程

void Gauss_jordan(){

/***** 行的交换&加减消元 *****/

for(re int i=1,r;i<=n;++i){ //正在处理第i行

r=i;

for(re int j=i+1;j<=n;++j)

if(fabs(a[j][i])>fabs(a[r][i])) r=j;

if(fabs(a[r][i])

puts("No Solution");return;

}

if(i!=r) swap(a[i],a[r]);

for(re int k=1;k<=n;++k){

//每一行都处理

if(k==i) continue;

double p=a[k][i]/a[i][i];

for(re int j=i;j<=n+1;++j) a[k][j]-=p*a[i][j];

}

}

//上述操作后会剩下对角矩阵,答案要除以系数

for(re int i=1;i<=n;++i) printf("%.2lf\n",a[i][n+1]/a[i][i]);

}

1.矩阵求逆

思路

求A的逆矩阵，把A和单位矩阵I放在一个矩阵里

对A进行加减消元使A化成单位矩阵

此时原来单位矩阵转化成逆矩阵

原理

A^{-1} * [AI] = [I A^{-1}]

举个栗子

求

\begin{bmatrix}2 & -1 & 0 \\-1 & 2 & -1 \\0 & -1 & 2\end{bmatrix}

首先

\begin{bmatrix}2 & -1 & 0 & 1 & 0 & 0 \\-1 & 2 & -1 & 0 & 1 & 0 \\0 & -1 & 2 & 0 & 0 & 1 \end{bmatrix}

对左边进行消元可得

\begin{bmatrix}2 & -1 & 0 & 1 & 0 & 0 \\ 0 & \frac{3}{2} & -1 & \frac{1}{2} & 1 & 0 \\0 & 0 & \frac{4}{3} & \frac{1}{3} &\frac{2}{3} & 1 \end{bmatrix}

此时已消成上三角矩阵，高斯消元开始回代，但约旦会消成对角矩阵

\begin{bmatrix}2 & 0 & 0 & \frac{3}{2} & 1 & \frac{1}{2} \\0 & \frac{3}{2} & 0 & \frac{3}{4} & \frac{3}{2} & \frac{3}{4} \\0 & 0 & \frac{4}{3} & \frac{1}{3} & \frac{2}{3} & 1\end{bmatrix}

最后每行除以系数

\begin{bmatrix}1 & 0 & 0 & \frac{3}{4} & \frac{1}{2} & \frac{1}{4} \\0 & 1 & 0 & \frac{1}{2} & 1 & \frac{1}{2} \\ 0 & 0 & 1 & \frac{1}{4} & \frac{1}{2} & \frac{3}{4} \end{bmatrix}

此时右半边即为所求

2.细节

开long long（不要冒风险，乘法很容易溢出）

模意义下除以一个数等于乘上逆元，可用快速幂求逆元（费马小定理）

Code

#include

#include

#include

#define re register

#define il inline

#define ll long long

using namespace std;

il ll read(){

ll s=0,f=0;char c=getchar();

while(c<'0'||c>'9') f=(c=='-'),c=getchar();

while(c>='0'&&c<='9') s=(s<<3)+(s<<1)+(c^'0'),c=getchar();

return f?-s:s;

}

const int N=405,mod=1e9+7;

int n;

ll a[N][N<<1];

il ll qpow(ll x,ll k){

ll ans=1;

while(k){

if(k&1) ans=ans*x%mod;

x=x*x%mod;

k>>=1;

}

return ans%mod;

}

il void Gauss_j(){

for(re int i=1,r;i<=n;++i){

r=i;

for(re int j=i+1;j<=n;++j)

if(a[j][i]>a[r][i]) r=j;

if(r!=i) swap(a[i],a[r]);

if(!a[i][i]){puts("No Solution");return;}

int kk=qpow(a[i][i],mod-2); //求逆元

for(re int k=1;k<=n;++k){

if(k==i) continue;

int p=a[k][i]*kk%mod;

for(re int j=i;j<=(n<<1);++j)

a[k][j]=((a[k][j]-p*a[i][j])%mod+mod)%mod;

}

for(re int j=1;j<=(n<<1);++j) a[i][j]=(a[i][j]*kk%mod);

//更新当前行 如果放在最后要再求一次逆元,不如直接放在这里

}

for(re int i=1;i<=n;++i){

for(re int j=n+1;j<(n<<1);++j) printf("%lld ",a[i][j]);

printf("%lld\n",a[i][n<<1]);

}

}

int main(){

n=read();

for(re int i=1;i<=n;++i)

for(re int j=1;j<=n;++j)

a[i][j]=read(),a[i][i+n]=1;

Gauss_j();

return 0;

}

网上浏览一圈头都要炸掉，线性代数太可怕了，定义好多

最后只看懂了这种方法

有什么问题欢迎评论区指出 ：）

参考文章

线性代数之——矩阵乘法和逆矩阵

逆矩阵的几种求法与解析(很全很经典)