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