学习线性代数心血来潮代码一波。
本文将简单介绍矩阵利用高斯消元完成的各种操作和科技,并用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.inv()
尚未定义。
1 |
|
OI选手可能会表示不满: 怎么都是double
啊?
取模请自己手动魔改…反正我现在不搞OI了,啦啦啦
显然,本文不会涉及辗转相除高斯消元。感兴趣者可以在OI资源分享中找到相关内容。
高斯消元,echelon form和reduced echelon form
初等行变换: 交换两行,一行乘等于一个常数,一行加等于若干倍的另一行。
高斯消元: 通过初等行变换将矩阵变成某个想要的形式。
高斯消元似乎也没啥好讲的,基本操作吧。大体思路就是模仿人手计算方式: 将除了第一个式子以外所有含有第一个变量的式子消掉第一个变量,将除了第二个式子以外所有含有第二个变量的式子消掉第二个变量…以此类推。
矩阵的(row) echelon form: 第行最左侧的非零元素的位置在第行最左侧的非零元素的位置的严格左侧。记表示第行最左侧的非零的下标,即(下标从0开始). 的(row) echelon form以后简记为.
矩阵的reduced (row) echelon form: 在 (row) echelon form的基础上满足. 的reduced (row) echelon form以后简记为.
高斯消元得到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; |
这里传入矩阵表示在将消成echelon form的同时对进行完全相同的操作。因为初等行变换等价于左乘一个矩阵,这里即分块矩阵的乘法分配: 。调用A.echelon(B)
,返回一个pair<Mat,Mat>
,即。
高斯消元得到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>
,即。
这里说明一下: double
高斯消元容易产生精度误差,一个优化的方法是每次选列主元的时候选择绝对值最大的。这里因为懒没有采用这种优化。
矩阵的迹、秩和逆
只有方阵才有定义。
矩阵的迹就是对角线上元素的和。这玩意有啥用?
消成echelon form以后非零行的数目就是矩阵的秩(rank)。
阶方阵可逆当且仅当。
. 对于一个可逆方阵来说,. 因此只要将捆绑着单位阵,消成,那么单位阵就会成为的逆。即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 _;} |
矩阵的行列式
只有方阵才有定义。
这里为了方便,暂时使用从开始的下标:
这里表示删除第行和第列后的矩阵。
这样不方便求,但是一行加等于另一行的变换不改变行列式、交换两行的变换使行列式取反,因此可以先高斯消元再求行列式。
而下或上三角矩阵行列式很好求,每次只考虑:
因此先转成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*_;} |
线性方程组
一个线性方程组可以表示为矩阵形式,其中和都是列向量。
有。而则是非常好解的: ,因此, 其余的()可以任取。
这里支持可视化输出线性方程组的通解: 给定已经转化为reduced echelon form的矩阵和,输出所有可能的。
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>
即,输出可能的的所有解(自由元独立输出,其余变量用自由元表示),具体形式可以自行试验,文末有试验。
LU分解和PLU分解
注意到每次解线性方程组复杂度都是高斯消元的复杂度。
如果每次更换,而不变,有没有更好的方法呢?
对于一个可逆方阵来说,即。可以直接求逆,求逆,此后每次计算只需要(矩乘列向量)。
然而对于非方阵怎么做呢?
U.P.D: 对于离线,可以直接构造,一次高斯消元解决。
LU分解是这样一种思路,设,则:
如果可以找到一一组合适的和,使得后两个方程很好解(),就做完了。
LU分解发现,为下三角矩阵,的时候似乎会更好解。
怎么分解呢?只需要高消就行了: ,而,整体即,因此可以捆绑高消求出然后直接求逆即可。
PLU分解几乎没有区别,只是为了方便: 如果分解过程中使用了交换操作,那么就会出现顺序问题,可以分解,其中是置换矩阵。
注: 并未代码实现过PLU分解,不做担保。
由于尚未实现PLU,因此在LU分解过程中必须保证不会发生交换。这是一个比较苛刻的条件。
留坑待填
考虑在做的时候发生了置换,那么这时候求出来的不再是,而是,其中是一个置换矩阵。
既然是一个下三角,那么通过选择排序就可以地转换为一个下三角。那么,就求出了和。
现在已经知道了,并且有,那么. 只需要每次把乘上(矩乘列向量当然也是的),然后再代入LU分解的做法即可。
将看做新的,此后只考虑LU分解。
如何计算?
如何计算?
这里如果要计算通解,就出现了一点问题,因为可能存在很多自由元。例如:
这时候首先可以得到. 如果任取自由元都代入当然没有问题,但是如果要求通解,保留整个公式的话,那么计算,这里的操作就不再是而是的了,因为本身就是一个无法化简的由多个自由元组成的表达式。
这样复杂度就再次退回到了(这里假设同阶)。
“GXZ分解"和"QGXZ分解”
首先不保证以下内容的正确性。其次不保证以下内容的独创性(希望大佬指出这东西叫啥)。
请教了隔壁宿舍大佬@GXZlegend。大佬给出了一种将LU分解的矩阵继续拆分来解决这个问题的方法。大致思路就是模仿高斯消元到底——反消。
设,其中是一个下三角矩阵(就是原来的),是一个上三角矩阵,。而其中就是原来的,就高斯消元的变换矩阵,事实上就是把高斯消元的正消和反消两步拆成了两个矩阵。
现在就有:
其中解和的方法和解原来相同,而现在,因此要输出通解——
当然是直接调用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大佬太强了!
首先解出,然后选择排序,于是就有:
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!!!
%%%%%%
扫描二维码即可在手机上查看这篇文章,或者转发二维码来分享这篇文章: