【线性代数】高斯消元与矩阵求逆

今天讲了线性代数，顺带复习了一下之前没有认真学的高斯消元以及矩阵求逆。

高斯消元：

考虑一个满秩的系数矩阵，它意味着有唯一解；而不存在唯一解的充要条件就是其行列式为 $0.$

那么考虑如何求解方程组：用初等行变换的形式将矩阵消成上三角矩阵，从而我们得到了最后一个未知数的解，再进行回代即可。

也可以直接消成一个对角矩阵。这个地方是和求逆一一样的。

用 $O(n^3)$ 可以完成高斯消元的过程，其中注意要确定用每次数绝对值最大的来消元，这样可以最大限度避免小数类计算的误差。

下面考虑矩阵求逆：

设已知矩阵为 $A,$ 目前要求 $P=A^{-1}.$

$$A\cdot A^{-1}=AP=I$$

那么如果我们通过 初等行变换 的形式，将 $A$ 消为 $I,$ 那么对左侧的 $I$ 进行同样的初等变换，对应地，它运算后的矩阵就应该是 $P$ 了。

于是同样在 $O(n^3)$ 的复杂度完成。

消元：

#include<bits/stdc++.h>
using namespace std;
const int N=105;
int n;
const double eps=1e-9;
double ans[N];
struct Matrix{
    double a[N][N];
    Matrix(){memset(a,0xfe,sizeof a);}
    void Guess(){
        for(int i=0;i<n;++i){
            int pos=i;
            for(int j=i+1;j<n;++j)
                if(fabs(a[pos][i])<fabs(a[j][i]))
                    pos=j;
            if(fabs(a[pos][i])<eps){
                puts("No Solution");
                exit(0);
            }
            if(i!=pos)swap(a[i],a[pos]);
            double div=a[i][i];
            for(int j=i;j<=n;++j)a[i][j]/=div;
            for(int j=i+1;j<n;++j){
                div=a[j][i];
                for(int k=i;k<=n;++k)a[j][k]-=div*a[i][k];
            }
        }
    }
}A;

int main(){
    scanf("%d",&n);
    for(int i=0;i<n;++i)
        for(int j=0;j<=n;++j)
            scanf("%lf",&A.a[i][j]);
    A.Guess();
    ans[n-1]=A.a[n-1][n];
    for(int i=n-2;~i;--i){
        ans[i]=A.a[i][n];
        for(int j=i+1;j<n;++j)ans[i]-=ans[j]*A.a[i][j];
    }
    for(int i=0;i<n;++i)printf("%.2lf\n",ans[i]);
    return 0;
}

求逆：

#include<bits/stdc++.h>
using namespace std;
const int N=500;
const int mod=1e9+7;
int n;
inline int Add(int x,int y){return (x+y)%mod;}
inline int Mul(int x,int y){return 1ll*x*y%mod;}
inline int Dec(int x,int y){return (x-y+mod)%mod;}
inline int Min(int x,int y){return x<y?x:y;}
inline int Max(int x,int y){return x>y?x:y;}
inline int Abs(int x){if(x<0)x=-x;return x;}
inline int qpow(int x,int y){
    int res=1;
    while(y){
        if(y&1)res=Mul(res,x);
        x=Mul(x,x);y>>=1;
    }
    return res;
}
inline int getinv(int x){return qpow(x,mod-2);}
struct Matrix{
    int a[N][N];
    Matrix(){memset(a,0,sizeof a);}
    void I(){for(int i=0;i<N;++i)a[i][i]=1;}
    void print(){
        for(int i=0;i<n;++i)
            for(int j=0;j<n;++j)
                printf("%d%c",a[i][j],j==(n-1)?'\n':' ');
    }
    void Init(){
        for(int i=0;i<n;++i)
            for(int j=0;j<n;++j)
                scanf("%d",&a[i][j]);
    }
    Matrix Getinv(){
        Matrix res;
        res.I();
        for(int i=0;i<n;++i){
            int pos=i;
            for(int j=i+1;j<n;++j)
                if(Abs(a[pos][i])<Abs(a[j][i]))
                    pos=j;
            if(!a[pos][i]){
                puts("No Solution");
                exit(0);
            }
            if(i!=pos)swap(a[i],a[pos]),swap(res.a[i],res.a[pos]);
            int div=a[i][i];
            div=getinv(div);
            for(int j=0;j<n;++j){
                res.a[i][j]=Mul(res.a[i][j],div);
                a[i][j]=Mul(a[i][j],div);
            }
            for(int j=0;j<n;++j){
                if(i==j)continue;
                div=a[j][i];
                for(int k=0;k<n;++k){
                    res.a[j][k]=Dec(res.a[j][k],Mul(div,res.a[i][k]));
                    a[j][k]=Dec(a[j][k],Mul(div,a[i][k]));
                }
            }
        }
        return res;
    }
}A;
int main(){
    scanf("%d",&n);
    A.Init();
    Matrix Ans=A.Getinv();
    Ans.print();
    // A.print();
    return 0;
}

关于应用：待补