【笔记】线性代数：矩阵的初等行变换

感觉这篇好乱哦……

关键词：高斯消元、初等行变换、行列式（的定义、性质、求法）、子式、代数余子式、伴随矩阵、克莱默法则、秩（的定义、性质、求法）、非奇异矩阵的逆、高斯消元的无解与无穷解区分。

目录

• 解方程 1：高斯消元
• 行列式
  • 行列式的定义
  • 行列式的性质
    • Lemma：排列奇偶性变换
    • Theorem 1：上/下三角矩阵
    • Theorem 2：初等行变换 1
    • Theorem 3：初等行变换 2
    • Theorem 4：行的线性性
    • Theorem 5：两行相同的矩阵
    • Theorem 6：初等行变换 3
  • 行列式的求值
  • 范德蒙德行列式
• 子式
  • \(k\) 阶子式
  • 代数余子式
  • 伴随矩阵
• 解方程 2：克莱默法则
• 秩
  • 秩的定义
  • 秩的求法
  • 秩的性质
  • 奇异矩阵
• 非奇异矩阵的逆
  • 伴随矩阵求法
  • 高斯消元求法
• 解方程 3：无解与无穷解判定
• 代码实现
  • 【模板】高斯消元
  • 【模板】矩阵求逆
  • 【模板】矩阵求逆（python）
  • 【模板】行列式求值

解方程 1：高斯消元

已知 \(n\) 元线性一次方程组。关于 \(x_1,x_2,\cdots,x_n\)。

\[\begin{cases} a_{1, 1} x_1 + a_{1, 2} x_2 + \cdots + a_{1, n} x_n = b_1 \\ a_{2, 1} x_1 + a_{2, 2} x_2 + \cdots + a_{2, n} x_n = b_2 \\ \cdots \\ a_{n,1} x_1 + a_{n, 2} x_2 + \cdots + a_{n, n} x_n = b_n \end{cases} \]

求出它的唯一解。或者报告无解、无数解。

考虑写成一个类似矩阵的形式（假设 \(n=3\)）：

\[\begin{bmatrix} a_{1,1} & a_{1,2} & a_{1,3} &|& b_1\\ a_{2,1} & a_{2,2} & a_{2,3} &|& b_2\\ a_{3,1} & a_{3,2} & a_{3,3} &|& b_3 \end{bmatrix} \]

对于这个矩阵我们可以有三种操作（即初等行变换），而不改变它的任何实质上的信息（即秩）：

• 交换任意两行。
• 对于一行的所有数字乘上同一个不为零的数。
• 将某一行的某个倍数与另一行对应位置相加，前者不变。

我们可以用这三种操作解这个矩阵，如果有解。我们的做法是这样的：

• 选一行的一个未知数。
• 化其为一。
• 用它消掉其他行的这个未知数，这样这个矩阵就只剩下它一个了。
• 重复上述操作 \(n\) 次以至于每一行都形如 \(x_i=b'_i\)。

没有唯一解的情况，就是找完所有行发现有一个未知数完全没有了。具体代码及实现细节参见最下方代码实现。

行列式

行列式的定义

行列式（Determinant）是对 \(n\) 阶方阵 \(A\) 定义的，是一个标量。\(A\) 的 \(n\) 阶行列式 \(\operatorname{det}(A)\) 或写作 \(|A|\) 定义如下：

\[\operatorname{det}(A)=\sum_p(-1)^{\operatorname{sgn}(p)}\prod_{i}A[i][p_i] \]

这里将排列的奇偶性定义为了 \(\operatorname{sgn}(p)\)，\(p\) 是排列，\(\operatorname{sgn}(p)\) 是 \(p\) 的逆序对个数的奇偶性。枚举所有排列，每一行挑一个乘起来，拿逆序对相关的东西做系数。

行列式的性质

Lemma：排列奇偶性变换

对于排列 \(p\)，交换 \(p_i,p_j\)（其中 \(1\leq i<j\leq n\)），\(\operatorname{sgn}(p)\) 的奇偶性改变。

假设 \(i<j\)，令 \(l=j-i\)。先将 \(a_j\) 调至 \(a_i\) 前，考虑这一步操作的影响。原来 \(a_j\) 与 \(a_{[i,j)}\) 的组合，顺序对设为 \(a\) 个，逆序对设为 \(b\) 个，那么 \(a+b=l\)，现在将 \(a_j\) 调至 \(a_i\) 前，顺序对变为 \(b\) 个，逆序对变为 \(a\) 个，所以变化量 \(a-b\equiv a-b+2b\equiv a+b\equiv l\pmod 2\)。再将 \(a_i\) 调至原来 \(a_j\) 的位置，变化量与 \(l-1\) 同奇偶，那么总的变化量是 \(l+l-1\equiv 1\pmod 2\)。

Theorem 1：上/下三角矩阵

单位矩阵的行列式为 \(1\)，上三角矩阵、下三角矩阵（对角线下/上全零的矩阵）的行列式为对角线乘积。

考虑选数乘起来的过程，如果不选零，那么归纳发现只有一种选法，就是选对角线。

Theorem 2：初等行变换 1

交换两行，行列式变号。

对于每个排列，这个交换两行只会影响 \(\operatorname{sgn}(p)\) 的奇偶性，由 Lemma 1 得 \(\operatorname{sgn}(p)\) 奇偶性改变，所以 \((-1)^{\operatorname{sgn}(p)}\) 全部变号，行列式变号。

Theorem 3：初等行变换 2

某一行乘以 \(t\)，行列式值乘以 \(t\)。

因为每行只选一个数进行乘积。

Theorem 4：行的线性性

由乘法分配律得到：

\[\begin{vmatrix}a_1+a_2&b_1+b_2\\c&d\end{vmatrix}=\begin{vmatrix}a_1&b_1\\c&d\end{vmatrix}+\begin{vmatrix}a_2&b_2\\c&d\end{vmatrix} \]

结合 Theorem 4 得到行的线性性，即这么一个东西：

\[\begin{vmatrix}a_1+\lambda a_2&b_1+\lambda b_2\\c&d\end{vmatrix}=\begin{vmatrix}a_1&b_1\\c&d\end{vmatrix}+\lambda \begin{vmatrix}a_2&b_2\\c&d\end{vmatrix} \]

意思是你可以对着行列式的一行线性地拆开。拆开这个动词是我现编的。

Theorem 5：两行相同的矩阵

有两行相同，行列式为 \(0\)。

如果交换这两行，行列式应该变号，但是方阵还是那个方阵，行列式的值应该没有改变。

Theorem 6：初等行变换 3

用一行的倍数加到另一行，行列式不变。

\[\begin{vmatrix}a&b\\c+\lambda a&d+\lambda b\end{vmatrix}=\begin{vmatrix}a&b\\c&d\end{vmatrix}+\lambda \begin{vmatrix}a&b\\a&b\end{vmatrix}=\begin{vmatrix}a&b\\c&d\end{vmatrix} \]

行列式的求值

Theorem 2,3,6 就是高斯消元用的三种操作，所以对方阵高斯消元，消成上三角矩阵，取其对角线乘积即可。

实现时，高斯消元要求有逆元；如果 \(\mathbb F_{mod}\) 下没有逆元，可以辗转消除消元。

范德蒙德行列式

这是一个行列式的例子。具体证明自己搜一下。

\[\begin{vmatrix} 1&1&\cdots&1\\ x_1&x_2&\cdots&x_n\\ x_1^2&x_2^2&\cdots&x_n^2\\ \vdots&\vdots&\ddots&\vdots\\ x_1^{n-1}&x_2^{n-1}&\cdots&x_n^{n-1}\\ \end{vmatrix}= \prod_{1\leq i<j\leq n}(x_j-x_i) \]

子式

\(k\) 阶子式

一个 \(n\times m\) 的矩阵 \(A\)，记 \(a_{i,j}\) 为 \(A\) 的第 \(i\) 行第 \(j\) 列的元素。定义记号 \(A_{S, T}\)，其中 \(S, T\) 是一个集合，\(|S|=|T|\)，首先求出一个 \(|S|\times |T|\) 的矩阵，由所有 \(i\in S, j\in T\) 的 \(a_{i, j}\) 组成，然后对其求行列式。

随机举个例子：\(A=\begin{bmatrix}1&2&3\\4&5&6\end{bmatrix}\)，取 \(S=\{1, 2\}, T=\{1, 3\}\)，则 \(A_{S, T}=\begin{vmatrix}1&3\\4&6\\\end{vmatrix}\)。

代数余子式

对于 \(n\times n\) 的方阵 \(A\)，定义 \(A\) 划去第 \(i\) 行第 \(j\) 列的行列式乘上 \((-1)^{i+j}\) 为 \(a_{i,j}\) 的代数余子式，记作 \(A_{i,j}=(-1)^{i+j}A_{[n]\setminus\{i\}, [n]\setminus\{j\}}\)。

有如下定理，不知道为什么：

1. \[\forall i, \sum_{j=1}^na_{ij}A_{ij}=|A| \]

2. \[\forall j, \sum_{i=1}^na_{ij}A_{ij}=|A| \]

3. \[\forall i\neq k, \sum_{j=1}^na_{ij}A_{kj}=0 \]

4. \[\forall j\neq k, \sum_{j=1}^na_{ij}A_{ik}=0 \]

伴随矩阵

对于 \(n\times n\) 的方阵 \(A\)，记 \(\operatorname{adj}(A)\) 为 \(A\) 的伴随矩阵，其第 \(i\) 行第 \(j\) 列元素为 \(A_{j, i}\)，即 \(a_{j, i}\) 的代数余子式。（注意伴随矩阵自带一个转置）

存在一种 \(O(n^3)\) 的计算伴随矩阵的算法，但是太难了，我不会，粘贴如下：

伴随矩阵求法

解方程 2：克莱默法则

又译作克莱姆法则。用于解形如 \(A\mathbf x=\mathbf b\) 的 \(n\) 元一次方程组，\(A\) 的大小为 \(n\times n\)，\(\mathbf x, \mathbf b\) 的大小都为 \(n\times 1\)，

克莱默法则：当 \(A\) 的行列式不为 \(0\) 时，\(x_i=\dfrac{|A_i|}{|A|}\)。其中 \(A_i\) 表示将 \(A\) 的第 \(i\) 列替换为 \(\mathbf b\) 得到的方阵，\(|A|\) 是 \(A\) 的行列式。

当 \(A\) 的行列式为 \(0\) 时，说明方程组无解或无穷解。如果已知方程组有至少一组解，同时 \(A\) 的行列式为 \(0\)，说明方程组无穷解。

据说另一种等价写法是 \(\mathbf x=A^{-1}\mathbf b\)。肯定是对的，但不知道为什么等价。

秩

秩的定义

秩有多种等价的定义。一种定义是说，一个 \(n\times m\) 的矩阵 \(A\) 的秩是某个 \(r\)，使得存在一个 \(r\) 阶子式非零，且若存在 \(r+1\) 阶子式则任意 \(r+1\) 阶子式都为零，这样的 \(r\) 称为矩阵 \(A\) 的秩，也记为 \(\operatorname{rank}(A)\) 或者好像有人写为 \(R(A)\)。

秩的求法

矩阵的初等行变换（交换两行、一行数乘非零数、一行的倍数加到另一行上）不改变矩阵的秩。

所以直接使用高斯消元求出矩阵的秩。

秩的性质

对于 \(n\times m\) 的矩阵 \(A\)，

• \(\operatorname{rank}(A)=\operatorname{rank}(A^T)\)，即其转置不改变其秩，也即行秩等于列秩。
• \(\operatorname{rank}(A)\leq \min(n, m)\)。
• \(\operatorname{rank}(\lambda A)=\operatorname{rank}(A)\ (\lambda\neq 0)\)。
• \(\operatorname{rank}(A)=0\iff A=0\)。
• \(\operatorname{rank}(AB)\leq \min(\operatorname{rank}(A), \operatorname{rank}(B))\)。
• \(\operatorname{rank}(A+B)\leq \operatorname{rank}(A)+\operatorname{rank}(B)\)。不知道为什么会这样。

奇异矩阵

对于 \(n\times n\) 的方阵 \(A\)，

• 若 \(R(A)=n\) 则称 \(A\) 为非奇异矩阵，推导出 \(|A|\neq 0\)，\(A\) 可逆。
• 若 \(R(A)<n\) 则称 \(A\) 为奇异矩阵，推导出 \(|A|=0\)，\(A\) 不可逆。

非奇异矩阵的逆

伴随矩阵求法

对于 \(n\times n\) 的方阵 \(A\)，其存在逆当且仅当 \(|A|\neq 0\)，其逆为 \(\frac{\operatorname{adj}(A)}{|A|}\)。

高斯消元求法

将一个单位矩阵 \(I_n\) 拼在 \(n\times n\) 方阵 \(A\) 的右边，并使用初等行变换将左半边消成单位矩阵，此时右半边就是 \(A\) 的逆。

为什么会这样？设所有初等行变换的矩阵的乘积为 \(P\)，则我们使得 \(PA=I_n\)，那么自然 \(P=A^{-1}\)，\(P\) 已经被刻画到右半边上了。

解方程 3：无解与无穷解判定

记 \(A|B\) 是把 \(A, B\) 这两个行数相等的矩阵横向拼在一起。求解 \(A\mathbf x=\mathbf b\) 时，

• 当 \(R(A)=R(A|\mathbf b)=n\)，原方程有唯一解。
• 当 \(R(A)=R(A|\mathbf b)<n\) 时，原方程有无穷解。
• 当 \(R(A)\neq R(A|\mathbf b)\) 时，原方程无解。（注意一定是 \(R(A)<R(A|\mathbf b)\) 了）

代码实现

【模板】高斯消元

P2455 [SDOI2006] 线性方程组 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)

Thanks https://www.luogu.com.cn/article/zje9tupn。

代码中使用 \(curr\) 维护当前已经做好了 \(curr\) 个元，现在正在尝试消第 \(curr + 1\) 个元。\(i\) 维护当前扫描到的列数，因为矩阵可能不满秩，所以 \(curr\) 可能不等于 \(i\)，需要分开维护。

每次循环中，首先看一下是否存在 \(curr\leq j< n\) 使得 \(a_{j, i}\neq 0\)。注意 \(j\geq curr\)，可能更小的 \(j\) 会有满足条件的，但是我们如果要了就会少掉之前消好的一个元，还不如不要。如果没有则使 \(i\) 自增 continue 掉，否则将这个 \(j\) 行交换到 \(curr\) 行来。接着将 \(a_{curr,i}\) 系数化为一，拿去对其他所有行做消元，并使 \(curr\) 自增即可。

使用 valarray 注意等号右边的东西不能在等号左边被改。意思是系数化为一时需要记一下原来是什么系数，否则可能会寄。

#include <bits/stdc++.h>
using namespace std;
#ifdef LOCAL
#define debug(...) fprintf(stderr, ##__VA_ARGS__)
#else
#define endl "\n"
#define debug(...) void(0)
#endif
using LL = long long;
constexpr double eps = 1e-3;
void simplify(valarray<double> a[], int n, int m) {
  int curr = 0;
  for (int i = 0; i < m && curr < n; i++) {
    int maxr = curr;
    for (int j = curr + 1; j < n; j++) {
      if (fabs(a[maxr][i]) < fabs(a[j][i])) maxr = j;
    }
    if (fabs(a[maxr][i]) < eps) continue;
    swap(a[maxr], a[curr]);
    double coe = a[curr][i];
    a[curr] /= coe;
    for (int j = 0; j < n; j++) if (j != curr) coe = a[j][i], a[j] -= a[curr] * coe;
    curr += 1;
  }
}
int getRank(valarray<double> a[], int n, int m) {
  int r = 0;
  for (int i = 0; i < n; i++) {
    bool flag = false;
    for (int j = 0; j < m; j++) flag |= fabs(a[i][j]) > eps;
    r += flag;
  }
  return r;
}
int n;
valarray<double> a[60];
int main() {
#ifndef LOCAL
  cin.tie(nullptr)->sync_with_stdio(false);  
#endif
  cin >> n;
  for (int i = 0; i < n; i++) {
    a[i].resize(n + 1);
    for (int j = 0; j <= n; j++) cin >> a[i][j];
  }
  simplify(a, n, n + 1);
  for (int i = 0; i < n; i++) {
    for (int j = 0; j <= n; j++) debug("%.2lf%c", a[i][j], " \n"[j == n]);
  }
  if (getRank(a, n, n) == n) {
    cout << fixed << setprecision(3);
    for (int i = 0; i < n; i++) cout << "x" << i + 1 << "=" << a[i][n] << endl;
  } else if (getRank(a, n, n + 1) == getRank(a, n, n)) {
    cout << 0 << endl; // 无穷解
  } else {
    cout << -1 << endl; // 无解
  }
  return 0;
}

【模板】矩阵求逆

P4783 【模板】矩阵求逆 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)

void simplify(valarray<mint> a[], int n, int m) {
  int curr = 0;
  for (int i = 0; i < m && curr < n; i++) {
    int maxr = curr;
    for (int j = curr + 1; j < n; j++) {
      if (a[j][i] != 0) maxr = j;
    }
    if (a[maxr][i] == 0) continue;
    swap(a[maxr], a[curr]);
    mint coe = a[curr][i];
    a[curr] /= coe;
    for (int j = 0; j < n; j++) if (j != curr) coe = a[j][i], a[j] -= a[curr] * coe;
    curr += 1;
  }
}
int getRank(valarray<mint> a[], int n, int m) {
  int r = 0;
  for (int i = 0; i < n; i++) {
    bool flag = false;
    for (int j = 0; j < m; j++) flag |= a[i][j] != 0;
    r += flag;
  }
  return r;
}
int n;
valarray<mint> a[410];
int main() {
#ifndef LOCAL
  cin.tie(nullptr)->sync_with_stdio(false);  
#endif
  cin >> n;
  for (int i = 0; i < n; i++) {
    a[i].resize(n * 2);
    for (int j = 0; j < n; j++) cin >> a[i][j].v;
    a[i][i + n] = 1;
  }
  simplify(a, n, n);
  for (int i = 0; i < n; i++) {
    for (int j = 0; j < n * 2; j++) debug("%d%c", a[i][j], " \n"[j + 1 == n * 2]);
  }
  if (getRank(a, n, n) == n) {
    for (int i = 0; i < n; i++) {
      for (int j = n; j < n * 2; j++) cout << a[i][j] << " \n"[j + 1 == n * 2];
    }
  } else {
    cout << "No Solution" << endl;
  }
  return 0;
}

【模板】矩阵求逆（python）

import numpy as np
matrix = np.array([
  [1, 1, 1, 1, 1, 1, 1, 1],
  [0, 1, 0, 1, 0, 1, 0, 1],
  [0, 0, 1, 1, 0, 0, 1, 1],
  [0, 0, 0, 1, 0, 0, 0, 1],
  [0, 0, 0, 0, 1, 1, 1, 1],
  [0, 0, 0, 0, 0, 1, 0, 1],
  [0, 0, 0, 0, 0, 0, 1, 1],
  [0, 0, 0, 0, 0, 0, 0, 1],
])
try:
    inverse_matrix = np.linalg.inv(matrix)
    print("矩阵的逆矩阵为：")
    print(inverse_matrix)
except np.linalg.LinAlgError:
    print("该矩阵不可逆。")

【模板】行列式求值

P7112 【模板】行列式求值 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)

const int P = 998244353;
LL mod(LL x) { return (x % P + P) % P; }
void red(LL& x) { x = mod(x); }
LL det(vector<LL>* a, int n) {
  LL res = 1;
  for (int i = 0; i < n; i++)
    for (LL& x : a[i]) x = mod(x);  //最好是正数
  for (int i = 0; i < n; i++) {
    for (int j = i + 1; j < n; j++) {
      while (a[i][i]) {  //最好是 a[i][i]，因为下面有个除法
        LL r = a[j][i] / a[i][i];  //整除
        for (int k = i; k < n; k++) red(a[j][k] -= a[i][k] * r);
        swap(a[i], a[j]), res = -res;
      }
      swap(a[i], a[j]), res = -res;
    }
  }
  for (int i = 0; i < n; i++) red(res *= a[i][i]);
  return mod(res);
}