学习线性代数心血来潮代码一波。
本文将简单介绍矩阵利用高斯消元完成的各种操作和科技，并用C++不使用任何第三方库手写一个封装有矩阵类型和多种操作的库。
主要内容包括: 矩阵的基本运算，矩阵的转置(transpose)，高斯消元，矩阵的echelon form和reduced echelon form，矩阵的迹(trace)，矩阵的秩(rank)，矩阵的逆(inverse)，矩阵的行列式(determinant)，高斯消元解线性方程组，线性方程组的通解，LU分解，PLU分解，以及由于我孤陋寡闻并没有查到的、大佬@GXZlegend传授的”GXZ分解”和”QGXZ分解”。
这里代码的实现更多采用了类似于”函数式编程”的模式，因此实现过程并没有很在意常数，常数可能比较大，使用时请谨慎。
反正我现在不搞OI了，暂时不用关心常数了，啦啦啦！
UPD 2019.10.25: 结尾处GXZ大佬本尊补充了另一种实现方式的”QGXZ分解”。更加OI，封装性相对较低，适合简便、快速地使用。
UPD 2019.10.27: 这种方法类似于LDU分解，两者形式上似乎不完全一样。
矩阵的基本运算和转置
加法、乘法、数乘、转置。似乎没啥好讲的。
下面定义基本的矩阵类型，封装了上述运算，以及矩阵的读入、输出、切割、右拼接。
特别地，这里给出了矩阵快速幂的重载运算符^，然而，为了通用性，快速幂的负数次方被定义为$A^{-K} = A.\text{inv()}^K$，目前这里的A.inv()尚未定义。
1 |
|
OI选手可能会表示不满: 怎么都是double啊？
取模请自己手动魔改……反正我现在不搞OI了，啦啦啦
显然，本文不会涉及辗转相除高斯消元。感兴趣者可以在OI资源分享中找到相关内容。
高斯消元，echelon form和reduced echelon form
初等行变换: 交换两行，一行乘等于一个常数，一行加等于若干倍的另一行。
高斯消元: 通过初等行变换将矩阵变成某个想要的形式。
高斯消元似乎也没啥好讲的，基本操作吧。大体思路就是模仿人手计算方式: 将除了第一个式子以外所有含有第一个变量的式子消掉第一个变量，将除了第二个式子以外所有含有第二个变量的式子消掉第二个变量……以此类推。
矩阵的(row) echelon form: 第$i$行最左侧的非零元素的位置在第$i+1$行最左侧的非零元素的位置的严格左侧。记$P(i)$表示第$i$行最左侧的非零的下标，即$P(i) \lt P(i+1), \forall i \in \lbrace 0,1,2,…,n-2 \rbrace$(下标从0开始). $A$的(row) echelon form以后简记为$A_{E}$.
矩阵的reduced (row) echelon form: 在 (row) echelon form的基础上满足$A{i,P(i)} = 1, A{k,P(i)} = 0, \forall i \in \lbrace 0,1,2,…,n-2 \rbrace, k \ne i$. $A$的reduced (row) echelon form以后简记为$A_{RE}$.
高斯消元得到echelon form(只正消不反消):
1 | inline pair<Mat,Mat> echelon(Mat B=Mat(0,0)){ assert(n==B.n||!B.n); Mat A(*this); A.d = 1; |
这里传入矩阵$B$表示在将$A$消成echelon form的同时对$B$进行完全相同的操作。因为初等行变换等价于左乘一个矩阵，这里即分块矩阵的乘法分配: $E(A \vert B) = (EA \vert EB)$。调用A.echelon(B)，返回一个pair<Mat,Mat>，即$(EA,EB)$。
高斯消元得到reduced echelon form(先正消后反消):
1 | inline pair<Mat,Mat> gauss(Mat B=Mat(0,0)){ assert(n==B.n||!B.n); Mat A(*this); |
同理，调用A.gauss(B)，返回一个pair<Mat,Mat>，即$(EA,EB)$。
这里说明一下: double高斯消元容易产生精度误差，一个优化的方法是每次选列主元的时候选择绝对值最大的。这里因为懒没有采用这种优化。
矩阵的迹、秩和逆
只有方阵才有定义。
矩阵的迹就是对角线上元素的和。这玩意有啥用？
消成echelon form以后非零行的数目就是矩阵的秩(rank)。
$n$阶方阵$A$可逆当且仅当$\text{rank}(A) = n$。
$A^{-1}(A \vert I) = (I \vert A^{-1})$. 对于一个可逆方阵来说，$I = A{RE}$. 因此只要将$A$捆绑着单位阵，消成$A{RE}$，那么单位阵就会成为$A$的逆。即A.inv()就是A.gauss(I).second。
1 | inline db trace(){assert(n==m); _ = 0; for(rint i = 0; i < n; _ += c[i][i], i++); return _;} |
矩阵的行列式
只有方阵才有定义。
这里为了方便，暂时使用从$1$开始的下标:
这里$\tilde A_{ij}$表示删除第$i$行和第$j$列后的矩阵。
这样不方便求，但是一行加等于另一行的变换不改变行列式、交换两行的变换使行列式取反，因此可以先高斯消元再求行列式。
而下或上三角矩阵行列式很好求，每次只考虑$i = 1$:
因此先转成echelon form，然后直接对角线乘积就好了。注意一下符号。
1 | inline db det(){assert(n==m); Mat A(this->echelon().X); _ = 1; for(rint i = 0; i < n; _ *= A(i,i), i++); return A.d*_;} |
线性方程组
一个线性方程组可以表示为矩阵形式$Ax = b$，其中$x$和$b$都是列向量。
有$E(A \vert b) = (A{RE} \vert b’)$。而$A{RE}x = b’$则是非常好解的: $A{i,P(i)} = 1$，因此$x{P(i)} = b’i - \sum{j \gt P(i)} x_j$, 其余的$x_k$($\forall i, k \ne P(i)$)可以任取。
这里支持可视化输出线性方程组的通解: 给定已经转化为reduced echelon form的矩阵$A_{RE}$和$b$，输出所有可能的$x$。
1 | inline void solution_set(pair<Mat,Mat> S){ int n = S.X.n, m = S.X.m; Mat A(S.X), b(S.Y); assert(n==b.n&&1==b.m); |
函数传入一个pair<Mat,Mat>即$(A_{RE},b)$，输出可能的$x$的所有解(自由元独立输出，其余变量用自由元表示)，具体形式可以自行试验，文末有试验。
LU分解和PLU分解
注意到每次解线性方程组复杂度都是高斯消元的复杂度$O(n^3)$。
如果每次$b$更换，而$A$不变，有没有更好的方法呢？
对于一个可逆方阵来说，$Ax = b$即$x = A^{-1}b$。可以直接求逆，求逆$O(n^3)$，此后每次计算$A^{-1}b$只需要$O(n^2)$(矩乘列向量)。
然而对于非方阵怎么做呢？
U.P.D: 对于离线，可以直接构造$(A \vert b_1 \vert b_2 … \vert b_q)$，一次高斯消元解决。
LU分解是这样一种思路，设$A = LU$，则:
如果可以找到一一组合适的$L$和$U$，使得后两个方程很好解($O(n^2)$)，就做完了。
LU分解发现，$L$为$n \times n$下三角矩阵，$U = A_E$的时候似乎会更好解。
怎么分解$A = LU = LA{E}$呢？只需要高消就行了: $E(A \vert I) = (A_E \vert E)$，而$A = LU = LA{E} = LEA$，整体即$L = E^{-1}$，因此可以捆绑高消求出$E$然后直接求逆即可。
PLU分解几乎没有区别，只是为了方便: 如果分解$A = LU$过程中使用了交换操作，那么就会出现顺序问题，可以分解$PA = LU$，其中$P$是置换矩阵。
注: 并未代码实现过PLU分解，不做担保。
由于尚未实现PLU，因此在LU分解过程中必须保证$A$不会发生交换。这是一个比较苛刻的条件。
留坑待填
考虑在做的时候发生了置换，那么这时候求出来的不再是$A = LU$，而是$A = (P^{-1}L)U$，其中$P$是一个置换矩阵。
既然$L$是一个下三角，那么$P^{-1}L$通过选择排序就可以$O(n^2)$地转换为一个下三角。那么$P(P^{-1}L \vert I) = (L \vert P)$，就求出了$P$和$L$。
现在已经知道了$P,L,U$，并且有$A = P^{-1}LU$，那么$Ax = b \Leftrightarrow LUx = Pb$. 只需要每次把$b$乘上$P$(矩乘列向量当然也是$O(n^2)$的)，然后再代入LU分解的做法即可。
将$Pb$看做新的$b$，此后只考虑LU分解。
如何计算$Ly = b$？
如何计算$Ux = y$？
这里如果要计算通解，就出现了一点问题，因为$U = A_E$可能存在很多自由元。例如:
这时候首先可以得到$x_4 = y_4 - 4(x_5+x_6+…+x_m)$. 如果任取自由元都代入$0$当然没有问题，但是如果要求通解，保留整个公式的话，那么计算$x_3 = y_3 - x_4 - 3(x_5+x_6+…+x_m)$，这里的$-x_4$操作就不再是$O(1)$而是$O(m)$的了，因为$x_4$本身就是一个无法化简的由多个自由元组成的表达式。
这样复杂度就再次退回到了$O(n^3)$(这里假设$n,m$同阶)。
“GXZ分解”和”QGXZ分解”
首先不保证以下内容的正确性。其次不保证以下内容的独创性(希望大佬指出这东西叫啥)。
请教了隔壁宿舍大佬@GXZlegend。大佬给出了一种将LU分解的$U$矩阵继续拆分来解决这个问题的方法。大致思路就是模仿高斯消元到底——反消。
设$A = GXZ$，其中$G$是一个$n \times n$下三角矩阵(就是原来的$L$)，$X$是一个$n \times n$上三角矩阵，$Z = A_{RE}$。而其中$XZ$就是原来的$U$，$GX$就高斯消元的变换矩阵$E$，事实上就是把高斯消元的正消和反消两步拆成了两个矩阵。
现在就有:
其中解$Gz = b$和$Xy = z$的方法和解原来$Ly = b$相同，而现在$Z = A_{RE}$，因此要输出通解——
当然是直接调用solution_set函数啊！
1 | inline void GXZ(Mat &G, Mat &X, Mat &Z){Mat I(n,n); for(rint i = 0; i < n; I(i,i) = 1, i++); pair<Mat,Mat> S(this->echelon(I)); G = S.Y.inv(); S = S.X.gauss(I); X = S.Y.inv(); Z = S.X;} |
1 | inline void GXZsolve(Mat G, Mat X, Mat Z, Mat b){ int n = G.n; assert(n==b.n&&1==b.m); Mat z(n,1), y(n,1), x(Z.m,1); db _; |
QGXZ分解其实就是PLU分解的思路应用到GXZ分解上。为什么换成了字母Q？，因为GXZ大佬太强了!
首先解出$A = (Q^{-1}G)XZ$，然后选择排序$Q(Q^{-1}G \vert I) = (G \vert Q)$，于是就有:
1 | inline void QGXZ(Mat &Q, Mat &G, Mat &X, Mat &Z){ |
1 | inline void QGXZsolve(Mat Q, Mat G, Mat X, Mat Z, Mat b){GXZsolve(G,X,Z,Q*b);} |
U.P.D: “QGXZ分解”的另一种实现
@GXZlegend大佬亲自实现的QGXZ分解。说明见[线性代数xOI/ACM]系数矩阵的QGXZ分解。
1 |
|
完整库以及测试文件
未经DEBUG慎用！！！
未经DEBUG慎用！！！
未经DEBUG慎用！！！
完整库:
1 |
|
测试文件:
1 |
|
输出结果:
1 | 1.000 1.000 0.000 0.000 |
总结
GXZx = Qb！！！
%%%%%%
扫描二维码即可在手机上查看这篇文章，或者转发二维码来分享这篇文章:
