高斯消元學習筆記
其實這個主題能夠復活主要還是粘了 \(\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[i][i]){
ll div=a[j][i]/a[i][i];
for(int k=i;k<=n;k++)a[j][k]=(a[j][k]-div*a[i][k])%P;
for(int k=i;k<=n;k++)swap(a[i][k],a[j][k]);
res=-res;
}
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;
}