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

一些无关紧要的Q&A

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

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

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

概述

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

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

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

前置技能

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

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

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

例如:

\[\begin{pmatrix} 1&1&0&0\\ 1&0&1&1\\ 2&1&0&0\end{pmatrix}= \begin{pmatrix} 1&0&0\\ 1&1&0\\ 2&0&1 \end{pmatrix} \begin{pmatrix} 1&1&0&0\\ 0&-1&1&1\\ 0&-1&0&0 \end{pmatrix}= \begin{pmatrix} 1&0&0\\ 1&1&0\\ 2&1&1 \end{pmatrix} \begin{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\times m\) 的矩阵,将其分解为一个 \(n\times n\) 的下三角矩阵 \(G\) 、一个 \(n\times n\) 的上三角矩阵 \(X\) 和一个 \(n\times m\) 的简化行阶梯矩阵(每个主元所在列的其它位置都是 \(0\) 的行阶梯矩阵) \(Z\) 的乘积的结果,即 \(A=G\times X\times Z\)

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

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

例如:

\[\begin{pmatrix} 1&1&0&0\\ 1&0&1&1\\ 2&1&0&0 \end{pmatrix}= \begin{pmatrix} 1&0&0\\ 1&1&0\\ 2&1&1 \end{pmatrix} \begin{pmatrix} 1&-1&0\\ 0&1&-1\\ 0&0&1 \end{pmatrix} \begin{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\times 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\) 行的前 \(\text{min}(2,3)-1\) 个数作为 \(L_1\) ,把 \(T_0\) 拿到 \(L_0\) 前面作为 \(L_1\) 即可。也即:

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

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

到此,解 \(Ax=b\) 就变为:分解 \(A=Q\times G\times X\times 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;
}
posted @ 2019-10-24 23:28  GXZlegend  阅读(906)  评论(4编辑  收藏