高斯消元学习笔记
高斯消元学习笔记
其实这个主题能够复活主要还是粘了 \(\text{LGV}\) 引理的光，不然我还不知道高斯消元其实不光能求解线性方程组。
求解线性方程组
这个只能说是典中典了，我不相信没有一个人的高斯消元不是从这里开始的。
我们考虑求解线性方程组的本质：将每一个式子所有未知数前都有系数转化成每一个式子只有一个未知数前有系数。那么我们不妨用矩阵来表示这样的关系。
对于线性方程组
我们将其写成
的形式，接着我们考虑消元，为了让整个矩阵最后呈现出
的形式，我们需要：
- 对于枚举到的每一行，将当前行对应的未知数前的系数变为 \(1\)。
- 将其他行对应此时的未知数前的系数变为 \(0\)。
最后从下往上，将每个未知数的值回代，求出每个未知数的解。
对于无解或无穷解的判断，主要在于全 \(0\) 行。如果出现 \(\begin{bmatrix}0&0&\cdots&0&0\end{bmatrix}\) 的行，那么说明有一个未知数可以取任意实数，因此有无穷解。如果出现 \(\begin{bmatrix}0&0&\cdots&0&b\end{bmatrix}\) 且 \(b\ne 0\)，那么说明 \(0\ne 0\)，那么此时整个方程组无解。复杂度是标准 \(O(n^3)\) 的。
代码
int gauss(){
int c,r;
for(c=1,r=1;c<=n;c++){
int t=r;
for(int i=r;i<=n;i++){
if(fabs(a[i][c])>fabs(a[t][c]))t=i;
}
if(fabs(a[t][c])<eps)continue;
for(int i=c;i<=n+1;i++)swap(a[t][i],a[r][i]);
for(int i=n+1;i>=c;i--)a[r][i]/=a[r][c];
for(int i=r+1;i<=n;i++){
if(fabs(a[i][c])>eps){
for(int j=n+1;j>=c;j--)a[i][j]-=a[r][j]*a[i][c];
}
}
r++;//消元
}
if(r<=n){
for(int i=r;i<=n;i++){
if(fabs(a[i][n+1])>eps)return 2;//方程无解
}
return 1;//方程有无穷多组解
}
for(int i=n;i>=1;i--){
for(int j=i+1;j<=n;j++)a[i][n+1]-=a[i][j]*a[j][n+1];
}
return 0;//方程有唯一解
}
求解矩阵逆元
学过矩阵快速幂的都知道，矩阵的单位元 \(I\) 形如：
参考整数逆元的定义：如果 \(aa^{-1}\equiv1\pmod p\)，则称 \(a^{-1}\) 是 \(a\) 在模 \(p\) 意义下的逆元，更一般的，如果没有模 \(p\)，那么 \(a^{-1}=\frac{1}{a}\)。
那么对于一个矩阵 \(A\)，是否存在它的逆元 \(A^{-1}\)，使得 \(AA^{-1}\equiv I\pmod p\) 或是 \(AA^{-1}=I\) 呢？答案是肯定的，但是矩阵的逆显然没有整数那么好求，于是我们考虑这样一个矩阵：
那么我们对左边的 \(A\) 进行消元，本质上就是将其变为单位元 \(I\)，也就是除以 \(A\)，那么矩阵将会变成 \(I|A^{-1}\)，此时右边的部分就是我们需要求的逆元了。
但是和整数逆元一样，矩阵在有些时候也没有逆元，比如不能成功除以 \(A\) 的时候，也就是消元失败的时候，当且仅当 \(\exists i,a_{i,i}=0\)，此时除数为 \(0\)，消元自然失败。因为和高斯消元的思路一致，复杂度是标准 \(O(n^3)\) 的。
我们可以考虑矩阵求逆和求解线性方程组的关系，不难发现，如果令
那么有 \(AX=B\)，即 \(X=A^{-1}B\)。由此我们发现，利用矩阵逆元我们也能求解出有解方程的正确解，但当矩阵没有逆元时，我们不能具体判断是无解还是有无穷多组解。
代码
int gauss(){
for(int i=1;i<=n;i++){
for(int j=i;j<=n;j++){
if(fabs(A[j][i])<=eps)continue;
for(int k=i;k<=(n<<1);k++)swap(A[i][k],A[j][k]);
break;
}
if(fabs(A[i][i])<=eps)return 1;
for(int j=(n<<1);j>=i;j--)A[i][j]/=A[i][i];
for(int j=1;j<=n;j++){
if(i==j)continue;
for(int k=(n<<1);k>=i;k--)A[j][k]-=A[j][i]*A[i][k];
}
}
return 0;
}
求解行列式
对于一个矩阵 \(A\) 的行列式，我们记为 \(\text{det}(A)\) 或 \(|A|\)。对于行列式有如下定义：
对于矩阵
有
其中 \(p\) 是 \(1\) 到 \(n\) 的排列，\(r(p)\) 是排列 \(p\) 中逆序对的个数。对于暴力搜索而言复杂度是 \(O(n!n\log n)\) 的。而对于行列式又有以下几个容易证得的性质：
- 如果矩阵 \(A\) 是一个上三角矩阵，即 \(\forall i>j,a_{i,j}=0\)，那么 \(|A|=\prod_{i=1}^{n}a_{i,i}\)。
- 如果矩阵 \(A\) 的某一行同时乘上一个常数 \(k\) 得到矩阵 \(A'\)，那么 \(|A'|=k|A|\)。
- 矩阵 \(A\) 的行列式与矩阵 \(A\) 的转置的行列式相等，即 \(|A|=|A^T|\)。
- 如果矩阵 \(A\) 的任意两行互换后得到矩阵 \(A'\)，那么 \(|A'|=-|A|\)。
- 如果矩阵 \(A\) 的其中一行加上另外一行得到矩阵 \(A'\)，那么 \(|A'|=|A|\)。
不难发现消元的过程符合性质 \(2,4,5\)，而消元后的矩阵符合性质 \(1\)，因此可以利用高斯消元求解。复杂度是 \(O(n^3)\) 的。
我们可以考虑行列式和求解线性方程组的关系，根据 \(\text{Cramar}\) 法则可以得知：
给定一个 \(n\) 元线性方程组，其未知数分别记作 \(x_1,x_2,\cdots,x_n\)。
设其矩阵表示形式为 \(AX=B\)，其中 \(A\) 为该方程组的增广矩阵，\(X\) 为由未知数组成的向量，\(B\) 为由常数项组成的向量，\(A_i\) 为 \(A\) 中的第 \(i\) 列被 \(B\) 取代后的矩阵，若 \(|A|\ne 0\)，则 \(x_i=\frac{|A_i|}{|A|}\)；否则该方程组无解。
由此，我们看出求解线性方程组也可以通过这样的方式完成。但是复杂度要更差。
代码
double gauss(){
double res=1;
for(int i=1;i<=n;i++){
for(int j=i;j<=n;j++){
if(fabs(A[j][i])<=eps)continue;
for(int k=i;k<=n;k++)swap(A[i][k],A[j][k]);
if(i!=j)res=-res;
break;
}
if(fabs(A[i][i])<=eps)return 0;
for(int j=i+1;j<=n;j++){
for(int k=n;k>=i;k--)A[j][k]-=A[j][i]/A[i][i]*A[i][k];
}
}
for(int i=1;i<=n;i++)res*=A[i][i];
return res;
}
//下面提供一种模数不为质数的写法
ll gauss(){
ll res=1;
for(int i=1;i<=n;i++){
for(int j=i+1;j<=n;j++){
while(a[j][i]){
ll div=a[i][i]/a[j][i];
for(int k=i;k<=n;k++)a[i][k]=(a[i][k]-div*a[j][k])%P;
for(int k=i;k<=n;k++)swap(a[i][k],a[j][k]);
res=-res;
}
}
}
for(int i=1;i<=n;i++)res=res*a[i][i]%P;
return res;
}
