[線性代數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;
}