一次线性方程组 高斯消元笔记

高斯消元原理

高斯消元用来解如下形式的方程组：

\[\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} \]

首先将所有系数写成系数矩阵：

\[\begin{bmatrix} a_{1, 1} & a_{1, 2} & \cdots & a_{1, n} & b_{1}\\ a_{2, 1} & a_{2, 2} & \cdots & a_{2, n} & b_{2}\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ a_{n, 1} & a_{n, 2} & \cdots & a_{n, n} & b_{n} \end{bmatrix} \]

接下来可以进行初等行列变换，将矩阵消成下三角矩阵。类似这样：

\[\begin{bmatrix} a'_{1, 1} & a'_{1, 2} & \cdots & a'_{1, n} & b'_{1}\\ 0 & a'_{2, 2} & \cdots & a'_{2, n} & b'_{2}\\ 0 & 0 & \cdots & a'_{3, n} & b'_{3}\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & \cdots & a'_{n, n} & b'_{n} \end{bmatrix} \]

然后最后一行就表示 \(a_{n, n} x_n = a_{n, n + 1}\)，可以解出 \(x_n\)。

然后往回带，解出 \(x_{n - 1}\)，以此类推即可。时间复杂度 \(O(n ^ 3)\)。

void gauss() {
	rep(i, 1, n) {
		rep(j, i, n) if (fabs(a[j][i]) > eps) {
			swap(a[j], a[i]); break;
		} if (fabs(a[i][i]) <= eps) continue;
		rep(j, i + 1, n) {
			if (fabs(a[j][i]) <= eps) continue;
			db inv = (db)a[j][i] / a[i][i];
			rep(k, i, n + 1)
				a[j][k] -= inv * a[i][k];
		}
	}
	for (int i = n; i; i -- ) {
		f[i] = a[i][n + 1] / a[i][i];
		for (int j = i - 1; j; j -- )
			a[j][n + 1] -= a[j][i] * f[i];
	}
}

高斯消元应用

• 解线性方程组

这个不用我说了吧。。。

• 在 \(\text{dp}\) 中用于消元

这个用处非常大。举个栗子：P3232 [HNOI2013] 游走

题意：给定 \(n\) 点 \(m\) 边无向简单图。要求给每条边重新编号 \(1 \sim m\)，使得每条边编号不同且从 \(1\) 到 \(n\) 经过路径的权值和期望尽可能小。

\(n \le 500, m \le 125000\)。

发现边数很多，所以从点数入手。设 \(f_i\) 表示第 \(i\) 个点被经过的期望次数，那么有：

\[\begin{cases} f_1 = \sum \limits_{(1, v) \in E, v \ne n}^{} \dfrac{f_v}{deg_v} + 1 \\ f_u = \sum \limits_{(u, v) \in E, v \ne n}^{} \dfrac{f_v}{deg_v} + 1 (u \ne 1)\\ \end{cases} \]

当然，\(n\) 号点被我们排除在外，因为到了 \(n\) 就停止了。

然后发现，上面的柿子形成了一个 \(n - 1\) 元方程。用高斯消元解方程就可以得到 \(f_u\)。

然后可以得到每条边被经过的期望次数：\(g_{(u, v)} = \dfrac{f_u}{deg_u} + \dfrac{f_v}{deg_v}\)。根据边被经过的期望次数贪心赋权即可。

其他用高斯消元求解 \(dp\) 的例题：P4321 随机漫游

• 矩阵求逆

矩阵 \(A\) 的逆矩阵定义为 \(A^{-1}\)，满足 \(A^{-1} \times A = I\) 且 \(A ^ {-1} \times A = I\)。其中 \(I\) 表示单位矩阵。

求法：将原矩阵 \(A\) 右面拼上一个单位矩阵 \(I\)。然后对左边一半也就是 \(A\) 做初等行列变换消成单位矩阵，同时对右边做左边同样的操作（左边行加右边也加，左边同乘右边也同乘）。最后得到的右半边就是 \(A ^ {-1}\)。

证明真的不会/kk

int gauss() {
	rep(i, 1, n) {
		rep(j, i, n)
			if (a[j][i] != 0) {
				swap(a[i], a[j]); break;
			}
		if (a[i][i] == 0) return 0;
		rep(j, i + 1, n) {
			if (a[j][i] == 0) continue;
			int inv = a[j][i] * qpow(a[i][i]) % mod;
			rep(k, i, 2 * n) a[j][k] -= inv * a[i][k] % mod,
				a[j][k] = (a[j][k] % mod + mod) % mod;
		}
	}
	rep(i, 1, n) {
		int inv = qpow(a[i][i]);
		rep(j, 1, 2 * n) a[i][j] = a[i][j] * inv % mod;
	}
	dep(i, n, 1) {
		rep(j, 1, i - 1) {
			int inv = a[j][i];
			rep(k, 1, 2 * n) a[j][k] -= a[i][k] * inv % mod,
				a[j][k] = (a[j][k] % mod + mod) % mod;
		}
	} return 1;
}
signed main() {
	read(n); int B = (10 * n + 5) * sizeof(LL);
	rep(i, 0, n + 1) a[i] = (LL*)malloc(B);
	rep(i, 1, n) rep(j, 1, n) read(a[i][j]);
	rep(i, 1, n) a[i][i + n] = 1; int s = gauss();
	if (!s) return puts("No Solution"), 0;
	for (int i = 1; i <= n; i ++ , pc('\n'))
		rep(j, 1, n) write(' ', a[i][j + n]);
	return 0;
}