[线性代数xOI/ACM]系数矩阵的QGXZ分解

一些无关紧要的Q&A

Q：你是怎么想到这个花里胡哨的算法的啊？
A：前几天学习线性代数时有幸和Magolor大佬讨论到 (LU) 分解在多解时的时间复杂度问题，于是yy出了这个奇怪（？）的算法。

Q：为什么叫 (QGXZ) 分解呀？你是不是在装逼啊？
A：这个名字是Magolor大佬起的，我也只能无条件服从咯~ 如有雷同绝非学术不端~

Q：Magolor大佬太强啦~
A：恭喜我们达成了共识~

概述

(QGXZ) 分解，是用于解决多线性方程组通解问题的算法。具体来讲：

给出 (n imes m) 的系数矩阵 (A) ，分别求 (Ax=b_1,Ax=b_2,...,Ax=b_q) 的通解 ，其中 (b_i) 是 (n imes 1) 的列向量。以下假设 (n,m,q) 同阶。

如果对 (b_i) 强制在线的话，朴素算法的时间复杂度为 (O(n^4)) 。如果对矩阵进行 (QGXZ) 分解，则复杂度降为 (O(n^3)) 。

前置技能

(QGXZ) 分解本质上是 (LU) 分解的扩展，因此先来介绍一下 (LU) 分解。

(LU) 分解是对于一个 (n imes m) 的矩阵，将其分解为一个 (n imes n) 的下三角矩阵 (L) 和一个 (n imes m) 的上梯形矩阵 (U) 的乘积的结果，即 (A=L imes U) 。

求法：对于矩阵 (A) ，从上到下进行矩阵行变换过程（这里仅考虑第三种行变换：将一行乘以一个数加到零一行上）。我们知道，使用一次行变换将 (A) 变成 (B) 的过程可以使用 (A=K imes B) 的形式描述，其中 (K) 是变换矩阵。由于在用上消下的前提下 (K) 是下三角矩阵，而下三角矩阵的乘积也是下三角矩阵，因此每次的变换矩阵的乘积就是我们所求的下三角矩阵 (L) ，而 (A) 的最终结果也是上梯形矩阵 (U)。

例如：

[egin{pmatrix} 1&1&0&0\ 1&0&1&1\ 2&1&0&0end{pmatrix}= egin{pmatrix} 1&0&0\ 1&1&0\ 2&0&1 end{pmatrix} egin{pmatrix} 1&1&0&0\ 0&-1&1&1\ 0&-1&0&0 end{pmatrix}= egin{pmatrix} 1&0&0\ 1&1&0\ 2&1&1 end{pmatrix} egin{pmatrix} 1&1&0&0\ 0&-1&1&1\ 0&0&-1&-1 end{pmatrix} ]

(LU) 分解有什么用？

假如现在有方程组 (Ax=b) ，它就等价于 (LUx=b) 。我们可以把 (Ux) 当作一个整体 (y) ，先解方程 (Ly=b) ，然后再解 (Ux=y) 。显然这两个方程都比较 “容易” 解出。

局限性

(LU) 分解有两点局限性：

1. 由于行变换的过程必须是使用上边的行消下边的行，因此对于一些矩阵可能不能直接进行 (LU) 分解；就算能进行 (LU) 分解，在处理小数时不能实现 “使用当前元系数绝对值最大的行消其余的行” ，精确度也就无法得到保证。

2. 即使矩阵能够进行 (LU) 分解，在解方程 (Ux=y) 时，如果方程有多解，则主元需要使用自由元来表示。而在代入求解的过程中，有 (O(n)) 个方程，每个方程要代入 (O(n)) 个主元，每个主元要用 (O(n)) 个自由元表示，因此就算知道了系数矩阵 (LU) 分解的形式，一次代入的复杂度也是 (O(n^3)) 的，和暴力没有区别。

下面我们介绍 (GXZ) 分解和 (QGXZ) 分解来解决这两点局限性。

(GXZ) 分解

(GXZ) 分解是对于一个 (n imes m) 的矩阵，将其分解为一个 (n imes n) 的下三角矩阵 (G) 、一个 (n imes n) 的上三角矩阵 (X) 和一个 (n imes m) 的简化行阶梯矩阵（每个主元所在列的其它位置都是 (0) 的行阶梯矩阵) (Z) 的乘积的结果，即 (A=G imes X imes Z) 。

这个求法也很简单：在LU分解使用行变换正消得到变换矩阵 (L) 和行阶梯矩阵 (U) 后，我们再反消一波，用主元行将上面行的相应位置消成 (0) ，并使用同 (LU) 分解的方法记录变换矩阵。由于每次都是用下面消上面，因此变换矩阵必然是上三角矩阵（和 (LU) 分解类似）。

在偷换一波变量名后便有 (A=GXZ) 。

例如：

[egin{pmatrix} 1&1&0&0\ 1&0&1&1\ 2&1&0&0 end{pmatrix}= egin{pmatrix} 1&0&0\ 1&1&0\ 2&1&1 end{pmatrix} egin{pmatrix} 1&-1&0\ 0&1&-1\ 0&0&1 end{pmatrix} egin{pmatrix} 1&0&0&0\ 0&-1&0&0\ 0&0&-1&-1 end{pmatrix} ]

这样的话，只需要解方程 (Gd=b) 、(Xe=d) 和 (Zx=e) 即可。前两个方程显然是 (O(n^2)) 的，而第三个方程只需要表示主元且没有代入过程，也是 (O(n^2)) 的。

于是我们就得到了一个 (O(n^3)) 预处理， (O(n^2)) 单次询问的算法。

(QGXZ) 分解

(GXZ) 分解处理了第二点局限性，第一点局限性则由 (QGXZ) 分解来解决。

(QGXZ) 分解即将 (n imes m) 的矩阵分解成置换矩阵 (Q) 和 (GXZ) 分解的乘积的形式。

具体方法：在 (GXZ) 分解的第一步（LU分解）时，假设当前已经消成了 (A=L_0U_0) 的形式，进一步变换消元时发现需要交换 (U_0) 的某两行，也即 (U_0=T_0U_1) ，其中 (T_0) 是置换矩阵。我们现在要做的就是将 (L_0T_0U_1) 变成 (T_1L_1U_1) ，即把 (L_0T_0) 变成 (T_1L_1) 。

我们知道，(L_0T_0) 相当于交换 (L_0) 的某两列，而 (T_1L_1) 相当于交换 (L_1) 的某两行。由于我们消元的过程是从上到下进行的，因此 (L_0) 要交换的两列必然是只有主对角线是 (1) ，其余位置为 (0) 。

因此，我们只需要手动交换 (L_0) 相应两行的主对角线前面的部分作为 (L_1) ，然后直接把 (T_0) 拿到前面，原封不动作为 (T_1) 即可。

例如：我们要交换 (L_0) 第 (2) 列和第 (3) 列，则手动交换 (L_0) 第 (2) 行和第 (3) 行的前 ( ext{min}(2,3)-1) 个数作为 (L_1) ，把 (T_0) 拿到 (L_0) 前面作为 (L_1) 即可。也即：

[egin{pmatrix} 1&0&0\ x&1&0\ y&0&1 end{pmatrix} egin{pmatrix} 1&0&0\ 0&0&1\ 0&1&0 end{pmatrix}= egin{pmatrix} 1&0&0\ 0&0&1\ 0&1&0 end{pmatrix} egin{pmatrix} 1&0&0\ y&1&0\ x&0&1 end{pmatrix} ]

每次交换都进行这样的过程，这样我们就把置换矩阵和置换矩阵放到了一起，把下三角矩阵和下三角矩阵放到了一起。由于它们的乘积都不会改变矩阵的特殊性质，因此最终的 (Q) 必然也是置换矩阵，(G) 必然也是下三角矩阵。

到此，解 (Ax=b) 就变为：分解 (A=Q imes G imes X imes Z) ，然后分别解 (Qc=b) 、(Gd=c) 、(Xe=d) 、(Zx=e) 即可。

单次询问的时间复杂度还是 (O(n^2)) 不变。

代码

老年选手不保证代码正确性（

#include <bits/stdc++.h>
#define N 510
#define eps 1e-6
using namespace std;
int pos[N];
double Q[N][N] , G[N][N] , X[N][N] , Z[N][N] , b[N] , c[N] , d[N] , e[N];
int main()
{
	int n , m , q , i , j , k , p = 0 , t;
	double mx;
	scanf("%d%d%d" , &n , &m , &q);
	for(i = 1 ; i <= n ; i ++ )
		for(j = 1 ; j <= m ; j ++ )
			scanf("%lf" , &Z[i][j]);
	for(i = 1 ; i <= n ; i ++ ) Q[i][i] = G[i][i] = X[i][i] = 1;
	for(i = 1 ; i <= m ; i ++ )
	{
		t = 0 , mx = eps;
		for(j = p + 1 ; j <= n ; j ++ )
			if(abs(Z[j][i]) > mx)
				t = j , mx = abs(Z[j][i]);
		if(!t) continue;
		pos[ ++ p] = i;
		for(k = i ; k <= m ; k ++ ) swap(Z[p][k] , Z[t][k]);
		for(k = 1 ; k <= n ; k ++ ) swap(Q[p][k] , Q[t][k]);
		for(k = 1 ; k < p ; k ++ ) swap(G[p][k] , G[t][k]);
		for(j = p + 1 ; j <= n ; j ++ )
		{
			G[j][p] = Z[j][i] / Z[p][i];
			for(k = i ; k <= m ; k ++ )
				Z[j][k] -= Z[p][k] * G[j][p];
		}
	}
	for(i = p ; i ; i -- )
	{
		for(j = i - 1 ; j ; j -- )
		{
			X[j][i] = Z[j][pos[i]] / Z[i][pos[i]];
			for(k = pos[i] ; k <= m ; k ++ )
				Z[j][k] -= Z[i][k] * X[j][i];
		}
	}
	while(q -- )
	{
		for(i = 1 ; i <= n ; i ++ ) scanf("%lf" , &b[i]);
		for(i = 1 ; i <= n ; i ++ )
			for(j = 1 ; j <= n ; j ++ )
				if(Q[i][j] == 1)
					c[j] = b[i];
		for(i = 1 ; i <= n ; i ++ )
		{
			d[i] = c[i];
			for(j = 1 ; j < i ; j ++ )
				d[i] -= G[i][j] * d[j];
		}
		for(i = n ; i ; i -- )
		{
			e[i] = d[i];
			for(j = n ; j > i ; j -- )
				e[i] -= X[i][j] * e[j];
		}
		for(i = p + 1 ; i <= n ; i ++ )
			if(abs(e[i]) > eps)
				break;
		if(i <= n) puts("No solution!");
		else
		{
			for(i = 1 ; i <= p ; i ++ )
			{
				printf("x[%d]=%lf" , pos[i] , e[i] / Z[i][pos[i]]);
				for(j = pos[i] + 1 ; j <= m ; j ++ )
					if(abs(Z[i][j]) > eps)
						printf("%+lfx[%d]" , -Z[i][j] / Z[i][pos[i]] , j);
				puts("");
			}
		}
	}
	return 0;
}