高斯消元学习笔记

标签：mat int 矩阵 long 学习笔记 using 对角高斯消

引入

高斯-约当消元法（Gauss–Jordan elimination）是求解线性方程组的经典算法，它在当代数学中有着重要的地位和价值，是线性代数课程教学的重要组成部分。

高斯消元法除了用于线性方程组求解外，还可以用于行列式计算、求矩阵的逆，以及其他计算机和工程方面。

过程

一个经典的问题，给定一个有 \(n\) 个未知数的线性方程组，求方程组的解。
例如：

\[\left\{ \begin{array}{} x_1 + 2x_2 = 4 \\ 2x_1 + 3x_2 = 5 \end{array} \right. \]

因为未知数的值在运算中也不会发生变化，可以尝试用矩阵将其转换。

定义

增广矩阵

一个矩阵，如果在矩阵右方有一个列向量或者一个列向量矩阵，那么这个矩阵就是增广矩阵。

将上方的例子改写成增广矩阵就是下图：

\[\left[ \begin{array} {c:|c} \begin{matrix} 1 & 2 \\ 2 & 3 \end{matrix} & \begin{matrix} 4\\ 5 \end{matrix} \end{array} \right] \]

可以发现最终的答案需要一个对角矩阵。

那么怎么将增广矩阵变换成对角矩阵呢。

将增广矩阵变成对角矩阵的过程

假定当前求解第 \(i\) 个未知数 \(x_i\)。

过程如下：

找到第一行 \(j\) 满足矩阵 \(mat\) 的 \(mat_{j,i} \neq 0 \land (mat_{j,j} \neq 1 \lor j \ge i)\) 。
交换第 \(i\) 行和第 \(j\) 行。
将其他行的第 \(i\) 列变成 \(0\)。
当 \(i \le n\) 时，重复以上步骤。

矩阵的秩

一个经过高斯消元的对角矩阵 \(A\) 如果满足第 \(i\) 行不是全零且第 \(i + 1\) 为全零，则矩阵的秩为 \(i\)，记作 \(r(A)\)。

无解情况

因为最后的结果会是一个对角方阵且此方阵的秩为 \(n\) 且满足此矩阵 \(A\) 中 \(\exists i, A_{i,i} = 0\) 并且右方常数不为 \(0\)。

无数解情况

满足一个对角矩阵 \(A\)，\(r(A) = n \land \exists i,A_{i,i} = 0\) 且其对应的常数为 \(0\)。

可以发现如果一个矩阵有唯一解满足 \(r(A) = n \land \forall i A_{i,i} \ne 0\)。

实现

一下代码为 SDOI 2006 线性方程组的代码。

#include <cstdio>
#include <cmath>
#include <algorithm>

using u32 = unsigned int ;
using i64 = long long ;
using u64 = unsigned long long ;
using f32 = float ;
using f64 = double ;

const f64 eps = 1e-6 ;
const int N = 105 ;

int n;
f64 mat[N][N];

void Gauss(){
    for (int i = 1; i <= n; i++) {
        int p = 0;
        for (int j = 1; j <= n; j++) {
            if (fabs(mat[j][j]) > eps && j < i) {
                continue;
            }

            if (fabs(mat[j][i]) > fabs(mat[p][i])) {
                p = j;
                break;
            }
        }

        if (!p) {
            continue;
        }
        std::swap(mat[i], mat[p]);

        for (int j = 1; j <= n; j++) {
            if (i == j) {
                continue;
            }

            double val = mat[j][i] / mat[i][i];
            for (int k = 1; k <= n + 1; k++) {
                mat[j][k] -= mat[i][k] * val;
            }
        }
    }
}

int main(){
    scanf("%d", &n);

    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n + 1; j++) {
            scanf("%lf", &mat[i][j]);
        }
    }

    Gauss();

    for (int i = 1; i <= n; i++) {
        if (fabs(mat[i][n + 1]) > eps && fabs(mat[i][i]) < eps) {
            puts("-1");
            return 0;
        }
    }
    
    for (int i = 1; i <= n; i++) {
        if (fabs(mat[i][n + 1]) < eps && fabs(mat[i][i]) < eps) {
            puts("0");
            return 0;
        }
    }

    for (int i = 1; i <= n; i++) {
        double ans = mat[i][n + 1] / mat[i][i];
        if (fabs(ans) < eps) {
            puts("0.00");
        } else {
            printf("%.2lf\n", ans);
        }
    }
    return 0;
}

时间复杂度

显然，\(\mathcal{O}(n ^ 3)\)。

其他应用

行列式求值

一个行列式如果它是对角矩阵，那么它的值为对角的乘积。

而高斯消元正好能将一个矩阵变换成对角矩阵，因此可以在 \(\mathcal{O}(n ^ 3)\) 求解。

实现

以下为【模板】行列式求值代码，通过辗转相减法实现了求值取模 \(p\) 的答案。

#include <cstdio>
#include <algorithm>

using i64 = long long ;

const int N = 605 ;

int n, mod;
i64 a[N][N];
i64 ans = 1;

void Gauss(){
    for (int i = 1; i <= n; i++) {
        for (int j = i + 1; j <= n; j++) {
            while (a[i][i]) {
                int t = a[j][i] / a[i][i] ;
                
                for (int k = 1; k <= n; k++) {
                    a[j][k] = (a[j][k] - t * a[i][k] % mod + mod) % mod;
                }

                ans *= -1;//行列式的性质，交换两行，需要改变结果的符号
                std::swap(a[i], a[j]);
            }

            ans *= - 1;
            std::swap(a[i], a[j]);
        }
    }
}

int main(){
    scanf("%d %d", &n, &mod);

    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            scanf("%lld", &a[i][j]);
        }
    }

    Gauss();

    for (int i = 1; i <= n; i++) {
        ans = (ans * a[i][i] + mod) % mod;
    }

    printf("%lld",(ans + mod) % mod);
    return 0;
}

矩阵求逆

矩阵求逆就是通过高斯消元将此矩阵消成单位矩阵，而右方可以记录矩阵的初等变换的积，最后右方就是矩阵的逆。

实现

以下为【模板】矩阵求逆的代码。

#include <cstdio>
#include <cstring>
#include <algorithm>

using u32 = unsigned int;
using i64 = long long ;
using u64 = unsigned long long ;

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

int n;
int mat[N][N << 1];

i64 qpow(i64 a, i64 b = mod - 2) {
    i64 ans = 1;
    
    for (; b; b >>= 1, a = a * a % mod) {
        if (b & 1) {
            ans = ans * a % mod;
        }
    }

    return ans;
}

int Gauss(){
    for (int i = 1; i <= n; i++) {
        int p = 0;
        for (int j = 1; j <= n; j++) {
            if (j < i && mat[j][j]) {
                continue;
            }
            if (mat[j][i]) {
                p = j;
                break;
            }
        }

        if (!p) {
            return -1;//如果方程组无解，自然就没有
        }
        std::swap(mat[i], mat[p]);

        i64 inv = qpow(mat[i][i]);
        for (int j = 1; j <= n; j++) {
            if (i == j) {
                continue;
            }

            i64 val = mat[j][i] * inv % mod;
            for (int k = 1; k <= (n << 1); k++) {
                mat[j][k] = (mat[j][k] - val * mat[i][k] % mod + mod) % mod;
            }
        }

        for (int j = 1; j <= (n << 1); j++) {
            mat[i][j] = (mat[i][j] * inv) % mod;
        }
    }

    return 1;
}

int main(){
    scanf("%d", &n);

    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            scanf("%d", &mat[i][j]);
            mat[i][n + i] = 1;
        }
    }

    if (Gauss() == -1) {
        printf("No Solution");
        return 0;
    }

    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            printf("%d ", mat[i][n + j]);
        }
        puts("");
    }
    return 0;
}

求解模意义下的方程组

很简单，将除法改成逆元就行了。

特殊的，二进制下可以用异或操作代替逆元。

标签：mat,int,矩阵,long,学习,笔记,using,对角,高斯消
From： https://www.cnblogs.com/zdrj/p/18224845

引入

过程

定义