时间戳

首页 > 技术文章 > 矩阵求逆·学习笔记 $\times$ [$LuoguP4783$]矩阵求逆

矩阵求逆·学习笔记 $\times$ [$LuoguP4783$]矩阵求逆

pks-t 2018-08-20 21:03 原文

哦？今天在\(luogu\)上fa♂现了矩阵求逆的板子……于是就切了切。

那么我们考虑一个矩阵\(A\)，它的逆矩阵记作\(A^{-1}\)，其中对于矩阵这个群来讲，会有\(A \cdot A^{-1} = I\) 其中\(I\)表示单位矩阵，主对角线均为\(1\) 。

那么我们对于矩阵\(A:\)

\(\begin{bmatrix} a_{1,1} & a_{1,2} & a_{1,3}\\ a_{2,1} & a_{2,2} & a_{2,3} \\ a_{3,1} & a_{3,2} & a_{3,3} \end{bmatrix} \quad\)

我们定义他的单位增广矩阵长底下这个样：

\(\begin{bmatrix} a_{1,1} & a_{1,2} & a_{1,3} & | & 1 & 0 & 0\\ a_{2,1} & a_{2,2} & a_{2,3} & | & 0 & 1 & 0 \\ a_{3,1} & a_{3,2} & a_{3,3} & | & 0 & 0 & 1 \end{bmatrix} \quad\)

最终这个矩阵应该是\(n\times 2m\)的。还有，中间那三条线是连起来的！连起来的！但是我并不知道怎么用\(\rm{LaTex}\)去搞这东西，于是\(GG\)

然后呢，我们对这个大矩阵进行暴力消元，直到发现左边的原矩阵变成\(I\)为止，右边的自然就是\(A^{-1}\) 那么我们思考它的正确性，以下是我想的根本不严谨只能意会的证明：

我们思考对于多个基于矩阵初等变换的变换，我们不妨称其为“变换集”。那我们思考对于原矩阵而言，从初态变成\(I\)的时候，经过的初等变换集虽然不唯一，但是一定组成一个变换集。那么如果右边的单位矩阵的变换集相同的话（意会证明来了），那么也就是说经历了左半个矩阵从原矩阵变成单位矩阵的全过程，那么原来的单位矩阵自然就变为了逆矩阵

好的，以上都是我不知道多长时间之前学的了\(233\)。啥？证明太水？没办法啊博主太弱啊……

那么思考代码：高消

\(hhh\)矩阵求逆·完

好吧，其实有个地方值得一说，就是我们朴素的高消都是消成上三角即可，包括求矩阵的秩，或者消元，或者求行列式。但是此处我们要求消成单位矩阵……嗯，这就很\(GG\).但是我从\(hjc\)那里听来一个很牛逼的思路：我们消到最后一定是可以把对角线消成\(1\)，这是很简单的，但，同时如果我们再用\(A_{i,i}\)去消同一列的其他元素，那么就一定可以将其消成\(0\)，因为此时窝萌的\(A_{i,i}\)是\(1\)啊就可以为所欲为了！那么也就是先正着消一遍，再反着消，用\(1\)元素去消其他的，然后就\(over\)了。那么大概效果图如下:

此处我们用\(F**k \ \ by \ \ A_{i,i}\)表示被\(1\)元素\(A_{i,i}\)消成\(0\)，那么

\(\begin{bmatrix} 1 & F**k \ \ by \ \ A_{2,4} & F**k \ \ by \ \ A_{3,3} & F**k \ \ by \ \ A_{4,4} \\ 0 & 1 & F**k \ \ by \ \ A_{3,3} &F**k \ \ by \ \ A_{4,4} \\ 0 & 0& 1 & F**k \ \ by \ \ A_{4,4} \\ 0 & 0 & 0 & 1\end{bmatrix}\)

\(233\)此处的\(F**k\)是一个你不可能见过的英文单词！不可能！！

上代码吧，此处的是\(LuoguP4783\)

#include <cstdio>
#include <iostream>
#define MAXN 401
#define ll long long
#define mod 1000000007

using namespace std ; 
ll NN, Inv ; bool flag ;
ll N, A[MAXN][MAXN << 1], i, j, k ;

inline ll qr(){
	ll k = 0, f = 1 ; char c = getchar() ;
	while (!isdigit(c)) {if (c == '-') f = -1 ; c = getchar() ;}
	while (isdigit(c)) k = (k << 3) + (k << 1) + c - 48, c = getchar() ;
	return k * f ; 
}
inline ll expow(ll P, ll Q){
	ll res = 1 ;
	while(Q){
		if (Q & 1) res = res * P % mod ; 
		P = P * P % mod ;
		Q >>= 1 ;
	}  return res ; 
}
inline void Gauss(){
	for (i = 1 ; i <= N ; ++ i){
		if (A[i][i] == 0) {flag = 1 ; return ;}
		Inv = expow(A[i][i], mod - 2) , A[i][i] = 1 ;
		for (j = i + 1 ; j <= NN ; ++ j) A[i][j] = A[i][j] * Inv % mod ;
		for (j = 1 ; j <= N ; ++ j)
			if (i == j) continue ;
			else {
				Inv = A[j][i], A[j][i] = 0 ;
				for (k = i + 1 ; k <= NN ; ++ k){
					A[j][k] = (A[j][k] - Inv * A[i][k]) % mod ;
					if (A[j][k] < 0) A[j][k] += mod ;
				}
			}
	}
}
int main(){
	cin >> N ; NN = N << 1 ; 
	for (i = 1; i <= N ; A[i][i + N] = 1, ++ i) 
		for (j = 1; j <= N ; ++ j) A[i][j] = qr() ;
	Gauss() ;
	if (flag) {
		cout << "No Solution" << endl ;
		return 0 ;
	}
	for (i = 1; i <= N ; ++ i){
		for(j = N + 1 ; j <= NN ; ++ j)
			printf("%lld ", A[i][j]) ;
		putchar('\n') ;
	}
	return 0 ;
}

\(Upd\)出锅了出锅了哈

上面的代码忘记\(swap\)了，所以很\(GG\)……

// luogu-judger-enable-o2
#include <cstdio>
#include <iostream>
#define MAXN 401
#define ll long long
#define mod 1000000007

using namespace std ; 
ll NN, Inv ; bool flag ;
ll N, A[MAXN][MAXN << 1], i, j, k ;

inline ll qr(){
    ll k = 0, f = 1 ; char c = getchar() ;
    while (!isdigit(c)) {if (c == '-') f = -1 ; c = getchar() ;}
    while (isdigit(c)) k = (k << 3) + (k << 1) + c - 48, c = getchar() ;
    return k * f ; 
}
inline ll expow(ll P, ll Q){
    ll res = 1 ;
    while(Q){
        if (Q & 1) res = res * P % mod ; 
        P = P * P % mod ;
        Q >>= 1 ;
    }  return res ; 
}
inline void Gauss(){
    for (i = 1 ; i <= N ; ++ i){
        for(j = i ; j <= N ; ++ j)
        if(A[j][i]){
            for(k = 1 ; k <= NN ; ++ k) swap(A[i][k], A[j][k]);  break ;
        }
        if (A[i][i] == 0) {flag = 1 ; return ;}
        Inv = expow(A[i][i], mod - 2) , A[i][i] = 1 ;
        for (j = i + 1 ; j <= NN ; ++ j) A[i][j] = A[i][j] * Inv % mod ;
        for (j = 1 ; j <= N ; ++ j)
            if (i == j) continue ;
            else {
                Inv = A[j][i], A[j][i] = 0 ;
                for (k = i + 1 ; k <= NN ; ++ k)
                    A[j][k] = ((A[j][k] - Inv * A[i][k]) % mod + mod) % mod ;
            }
    }
}
int main(){
    cin >> N ; NN = N << 1 ;
    for (i = 1; i <= N ; ++ i) 
        for (j = 1; j <= N ; ++ j)
            A[i][j] = qr() ;
    for (i = 1; i <= N ; ++ i) A[i][i + N] = 1 ;
    Gauss() ;
    if (flag) {
        cout << "No Solution" << endl ;
        return 0 ;
    }
    for (i = 1; i <= N ; ++ i){
        for(j = N + 1 ; j <= NN ; ++ j)
            printf("%lld ", A[i][j]) ;
        putchar('\n') ;
    }
}

推荐阅读