Linear Equation - 演算法筆記
文章推薦指數： 80 %
由上往下消除變數的過程，就叫做「高斯消去法」。
LinearEquation
LinearEquation
「一次方程式」。
一次函數的等式。
3x+4y-6=2x+5z-1
移項整理：變數在左邊，常數在右邊。
1x+2y-5z=5
SystemofLinearEquations
「一次方程組」。
許多道一次方程式同時成立。
{1x+2y-5z=5
{2x+4y+6z=1
{3x+1y+7z=4
一次方程組求解，等同於線性函數求解，或者說矩陣求解。
{1x+2y-5z=5[12-5][x][5]
{2x+4y+6z=1--->[246][y]=[1]
{3x+1y+7z=4[317][z][4]
Axb
特殊矩陣擁有特殊演算法，稍微省時，詳情請見MATLAB的mldivide說明圖片。
簡單起見，以下不介紹特殊矩陣的演算法。
Inverse
「反矩陣」。
一旦知道反矩陣，即可輕鬆解一次方程組。
而且還可以隨意抽換等號右邊數值。
solveAx=b--->findA⁻¹,thenx=A⁻¹b
一次方程組求解演算法，稍作修改，即得反矩陣演算法。
LinearEquation:Elimination
等量公理
相信大家已經學過如何解一次方程組：由上往下消除變數，變成階梯狀；由下往上解出變數，變成對角線。
由上往下消除變數的過程，就叫做「高斯消去法」。
演算法（GaussianElimination）
一個矩陣，化成上三角矩陣，對角線元素皆為一。
由上往下，處理每個橫條，實施三種運算：交換、倍率、相減。
如果橫條的首項係數不是零，以該橫條抵銷下方橫條，令下方橫條的首項係數化成零。
如果是零，就往下找到非零橫條，先交換，再抵銷。
pivotrow：首項係數不是零的橫條。
pivot：上述橫條的首項係數。
pivoting：找到pivotrow，視情況進行交換。
減少誤差的方法：取絕對值最大的pivotrow，抵銷其餘row。
換句話說，首項係數絕對值最大的橫條，總是交換到上方，再來抵銷下方橫條。
倍率小於1，相減誤差少。
矩陣邊長N×M，時間複雜度O(N²M)。
大家習慣討論方陣，N=M，時間複雜度O(N³)。
方陣的高斯消去法，程式碼如下所示。
矩陣的高斯消去法，留給大家自行練習。
doublematrix[10][10];
voidgaussian_elimination()
{
for(inti=0;i<10;++i)
{
//如果橫條的首項係數為零，嘗試與下方橫條交換。
if(matrix[i][i]==0)
for(intj=i+1;j<10;++j)
if(matrix[j][i]!=0)
{
//交換上方橫條與下方橫條。
for(intk=i;k<10;++k)
swap(matrix[i][k],matrix[j][k]);
break;
}
//如果首項係數都是零，那就略過。
if(matrix[i][i]==0)continue;
//橫條的首項係數調整成一。
為了讓對角線皆為一。
doublet=matrix[i][i]; //首項係數
for(intk=i;k<10;++k)
matrix[i][k]/=t;
//抵銷下方橫條，令下方橫條的首項係數化成零。
for(intj=i+1;j<10;++j)
if(matrix[j][i]!=0)
{
doublet=matrix[j][i];
for(intk=i;k<10;++k)
matrix[j][k]-=matrix[i][k]*t;
}
}
}
解一次方程組
首先實施高斯消去法，求得上三角矩陣。
由下往上，處理每個橫條。
把先前解出的變數，代入到目前橫條，解出變數。
時間複雜度O(N²)。
doublematrix[10][10+1]; //參數化的矩陣
doublex[10]; //一次方程組的解
voidback_substitution()
{
for(inti=10-1;i>=0;--i)
{
if(matrix[i][i]==0)return; //無解、多解
doublet=0;
for(intk=i+1;k<10;++k)
t+=matrix[i][k]*x[k];
x[i]=(matrix[i][10]-t)/matrix[i][i];
}
}
UVa101091052410828ICPC3563
演算法（Gauss–JordanElimination）
「高斯喬登消去法」是延伸版本。
對角線化成一、其餘化成零。
時間複雜度仍是O(N³)。
高斯喬登消去法也可以解一次方程組，但是步驟數量比較多。
主要用途是求反矩陣。
求反矩陣
利用「高斯喬登消去法」。
doublematrix[10][10*2]; //參數化的矩陣
boolgauss_jordan_elimination()
{
//填好參數化的部分
for(inti=0;i<10;++i)
for(intj=0;j<10;++j)
matrix[i][10+j]=0;
for(inti=0;i<10;++i)
matrix[i][10+i]=1;
//開始進行高斯喬登消去法
//內容幾乎與高斯消去法相同
for(inti=0;i<10;++i)
{
if(matrix[i][i]==0)
for(intj=i+1;j<10;++j)
if(matrix[j][i]!=0)
{
for(intk=i;k<10*2;++k)
swap(matrix[i][k],matrix[j][k]);
break;
}
//反矩陣不存在。
if(matrix[i][i]==0)returnfalse;
doublet=matrix[i][i];
for(intk=i;k<10*2;++k)
matrix[i][k]/=t;
//消去時，所有橫條通通消去。
for(intj=0;j<10;++j)
if(i!=j&&matrix[j][i]!=0)
{
doublet=matrix[j][i];
for(intk=i;k<10*2;++k)
matrix[j][k]-=matrix[i][k]*t;
}
}
returntrue;
}
求determinant
利用「高斯消去法」，對角線不化成一，保留原有數字。
上三角矩陣，對角線元素的乘積，便是determinant。
如果矩陣裡都是整數，那麼determinant也是整數。
想避免浮點數誤差，可以使用輾轉相除法進行消去。
時間複雜度O(N³logC)，C是絕對值最大的首項係數。
intmatrix[10][10];
voidrowgcd(int*a,int*b,inti)
{
if(abs(a[i])
高斯消去法直接得到正解。
鬆弛法逐步逼近正解。
當b=0，則完全等價於D⁻¹(L+U)實施不動點遞推法。
判斷收斂：檢查D⁻¹(L+U)的特徵值的絕對值是否都小於1。
「嚴格對角優勢strictlydiagonallydominant」：每個橫條，對角線元素大於其餘元素總和；每個直條亦然。
滿足sdd，保證收斂；不滿足時，可能收斂、也可能不收斂。
foreachrow,|Aii|>∑|Aij|
j≠i
演算法（Gauss–SeidelIteration）
每回合依序計算xyz，剛出爐的數字，馬上拿來使用，加快收斂速度。
[43-1][x][1]
[251][y]=[2]
[-2-26][z][3]
{4x+3y-z=1{x=(1-3y+z)/4
{2x+5y+z=2=>{y=(2-2x-z)/5
{-2x-2y+6z=3{z=(3+2x+2y)/6
[x₀][0]
[y₀]=[0]隨便設定一個初始值
[z₀][0]
[x₁][(1-3y₀+z₀)/4]依序計算x₁、y₁、z₁，
[y₁]=[(2-2x₁-z₀)/5]剛出爐的數字，馬上拿來使用，
[z₁][(3+2x₁+2y₁)/6]加快收斂速度。
xₖ₊₁=D⁻¹(b-Uxₖ-Lxₖ₊₁)
演算法（SuccessiveOver-relaxation）
原數值、新數值，以固定比例混合。
可令原數值權重小於零、新數值權重大於一，硬是催出收斂速度。
[x₁][(1-w)⋅x₀+w⋅(1-3y₀+z₀)/4]
[y₁]=[(1-w)⋅y₀+w⋅(2-2x₁-z₀)/5]
[z₁][(1-w)⋅z₀+w⋅(3+2x₁+2y₁)/6]
LinearEquation:Projection
LinearEquation與Geometry
一次方程式得視作幾何元件：點、線、面、……。
ContourPlot3D[1x+2y-5z==5,{x,-5,5},{y,-5,5},{z,-5,5},Boxed->False,Axes->False,Mesh->None]
一次方程組的解得視作一堆幾何元件的交集。
f:=1x+2y-5z-5;g:=2x+4y+6z-1;h:=3x+1y+7z-4;ContourPlot3D[{f==0,g==0,h==0},{x,-5,5},{y,-5,5},{z,-5,5},Boxed->False,Axes->False,Mesh->None,ContourStyle->Directive[Opacity[0.5]]]
演算法（Kaczmarz'sMethod）
隨便設定一解，依序投影到第一道等式、第二道等式、……，不斷循環，逐步逼近正解。
時間複雜度O(N²T)，N是方陣維度，T是遞推次數。
無解時，最後形成無窮迴圈，稱作「極限環LimitCycle」。
https://books.google.com.tw/books?id=1PJ-WHepeBsC&pg=PA620
數學公式（Cramer'sRule）
兩線交點的演算法：求得平行四邊形的面積，以面積比例求得交點位置。
請見本站文件「Intersection」。
「克拉瑪公式」則是此演算法的高維度版本，形成了非常漂亮的公式解！求得超平行體的容積，以容積比例求得解。
linearequation:
Ax=b
solution:
{x=det(Aˣ)/det(A)
{y=det(Aʸ)/det(A)
{z=det(Aᶻ)/det(A)
unisolvence:
Ax=bhasauniquesolutioniffdet(A)≠0
example:
{1x+2y-5z=5[12-5][x][5]
{2x+4y+6z=1--->[246][y]=[1]
{3x+1y+7z=4[317][z][4]
Axb
_b_b_b
[|5|2-5][1|5|-5][12|5|]
Aˣ=[|1|46]Aʸ=[2|1|6]Aᶻ=[24|1|]
[|4|17][3|4|7][31|4|]
‾‾‾
determinant是矩陣當中所有向量所構成的超平行體的容積。
時間複雜度等於N+1次determinant的時間複雜度，O(N⁴)。
Determinant
determinant起初用來判定一個一次方程組是否有解、解是多少，因而稱作「決定因子」。
古人沒有意識到determinant是容積。
字面意義是「決定因子」，中文教科書卻譯作「行列式」。
真是異想天開的翻譯啊！
http://mathworld.wolfram.com/DeterminantExpansionbyMinors.html
行列式的計算過程是：先刪除一橫行，接著分別刪除每一直行，形成N-1個(N-1)×(N-1)子矩陣，添上正負號。
原矩陣的行列式，等於這些子矩陣的行列式總和。
每個子矩陣各自遞迴下去，直到N=1。
1×1矩陣的行列式，等於矩陣元素。
時間複雜度O(N!)。
N=2or3的時候比較特別，可以直接累加所有「左上右下斜線」的乘積、累減「右上左下斜線」的乘積。
中學數學課程有教。
計算行列式，也可以使用高斯消去法，時間複雜度O(N³)。
不過與其採用高斯消去法求行列式、再用行列式解一次方程組，不如直接採用高斯消去法解一次方程組。
就當作是學個想法吧。
LinearEquation:Decomposition
LinearEquation與Algebra
一次方程組可以寫成矩陣乘法的形式。
{1x+2y-5z=5[12-5][x][5]
{2x+4y+6z=1--->[246][y]=[1]
{3x+1y+7z=4[317][z][4]
Axb
一次方程組得視作線性函數，解一次方程組得視作線性反函數。
-1
[12-5][x][5][x][12-5][5]
[246][y]=[1]--->[y]=[246][1]
[317][z][4][z][317][4]
AxbxA⁻¹b
數學公式（InverseMatrix）
一旦求得反矩陣，即可輕鬆解一次方程組。
solveAx=b--->findA⁻¹,thenx=A⁻¹b
計算反矩陣，使用高斯喬登消去法，時間複雜度O(N³)。
不過與其採用高斯喬登消去法求反矩陣、再用反矩陣解一次方程組，不如直接採用高斯消去法解一次方程組。
就當作是學個想法吧。
數學公式（Eigendecomposition）
一旦求得特徵分解，即可輕鬆解一次方程組。
A=EΛE⁻¹
A⁻¹=EΛ⁻¹E⁻¹
x=A⁻¹b=EΛ⁻¹E⁻¹b
不過與其採用特徵分解求反矩陣、再用反矩陣解一次方程組，不如直接採用高斯消去法解一次方程組。
就當作是學個想法吧。
LinearEquation:RootFinding
演算法（Newton'sMethod）
Ax=b的解，即是f(x)=Ax-b的根。
牛頓法。
x=x-Aᵀ⁻¹(Ax-b)。
光是求反矩陣就需時O(N³)。
無人使用。
不過與其求反矩陣、不斷遞推，不如直接採用高斯消去法。
就當作是學個想法吧。
演算法（RichardsonIteration）
Ax=b的解，即是g(x)=Ax-b+x的不動點。
不動點遞推法。
可以添加權重α。
x=x+α(Ax-b)。
時間複雜度O(N²T)，N是方陣維度，T是遞推次數。
不過與其求餘數、不斷遞推，不如直接採用JacobiIteration。
就當作是學個想法吧。
LinearEquation:Optimization
「一次函數求解」等於「二次函數求極值」
A必須是對稱正定矩陣，使得二次函數有唯一最小值。
minimizef(x)=1/2xᵀAx-bᵀx
solvef′(x)=Ax-b=0
solveAx=b
f(x)的最小值位置，即是f′(x)的根，即是Ax=b的解。
f(x)的最佳化演算法是「共軛梯度法」。
優於高斯消去法。
時間複雜度O(N²T)，N是方陣維度，T是遞推次數。
「矩陣求解」化作「對稱正定矩陣求解」
當A不是對稱正定矩陣，那麼等號兩邊同時乘上Aᵀ，得到對稱半正定矩陣AᵀA。
當A有唯一解，則得到對稱正定矩陣AᵀA。
solveAx=b(ifAhasuniquesolution)
--->solveAᵀAx=Aᵀb(AᵀAissymmetricpositivedefinite)
--->solveA'x=b'(uniquesolution)
solveAx=b(ifAhasnosolution)
--->min‖Ax-b‖²(leastsquares)
--->solveAᵀAx=Aᵀb(normalequation)
--->solveA'x=b'(leastsquaressolution)
不過與其變成對稱正定矩陣、再用共軛梯度法解一次方程組，不如直接採用高斯消去法解一次方程組。
就當作是學個想法吧。
LinearLeastSquares
LinearLeastSquares
唯一解是稀奇的，無解、多解是普遍的。
無解、多解時，改為找到平方誤差最小的解。
nosolution：無解。
改求‖Ax-b‖²最小的解。
infinitelymanysolutions：多解。
改求‖x‖²最小的解。
uniquesolution：唯一解。
求Ax=b的解。
min‖Ax-b‖²：方程組每一道等式，求得等號左右兩邊的差的平方；累計所有等式，總和越小越好。
min‖x‖²：解的每一項的平方，總和越小越好。
Least意指「盡量小」，Squares意指「平方和」。
平方誤差比起絕對值誤差，有三個好處：一、讓每道等式的誤差保持均勻，不會有某道等式誤差特別高。
二、「一次微分等於零」容易推導數學公式。
三、誤差函數是開口向上的拋物線函數，有唯一最小值，沒有鞍點，容易最佳化。
平方誤差，誤差高低可以拿來比較，適合成為一個參考指標。
NoSolution/InfinitelyManySolutions
請見「Rouché-CapelliTheorem」。
矩陣的rank、column，得以區分無解、多解。
rank：消去冗餘的等式。
column：直條數量、變數數量。
nosolution:rank(A)≠rank([A|b])
infinitelymanysolutions:rank(A)=rank([A|b])≠column(A)
uniquesolution:rank(A)=rank([A|b])=column(A)
OverdeterminedSystem/UnderdeterminedSystem
[XX][XXX]
[XX][XXX][XXX]
[XX][XXX][XXX]
thinfatsquare
overdeterminedsystem(thinmatrix):row(A)>column(A)
underdeterminedsystem(fatmatrix):row(A)
「一次微分等於零」的地方是極值、鞍點。
因為平方誤差是開口向上的拋物面，所以「一次微分等於零」的地方必是最小值。
以下只證明無解的情況。
時間複雜度O(N³)。
solveAx=b
min‖Ax-b‖²
d/dx‖Ax-b‖²=0「一次微分等於零」的地方是極值、鞍點
[d/dx(Ax-b)][2(Ax-b)]=0微分連鎖律
Aᵀ[2(Ax-b)]=0微分
AᵀAx=Aᵀb同除以2、展開、移項
x=(AᵀA)⁻¹Aᵀb移項。
注意到A的向量們必須線性獨立！
注意到最後一步，A的向量們必須線性獨立（事先清除冗餘的、無意義的變數），AᵀA才有反矩陣。
結果宛如向量投影公式：b投影到A。
Aᵀbdot(A,b)
x=(AᵀA)⁻¹Aᵀb=―――=――――――――=projAb
AᵀAdot(A,A)
平方誤差盡量小＝垂直投影！分析與幾何兩大領域打通了！
solveAx=b
--->min‖Ax-b‖²
--->2Aᵀ(Ax-b)=0
--->x=projAb
數學公式（QRDecompostion）
http://www.cs.cornell.edu/courses/cs322/2007sp/notes/qr.pdf
以下只證明無解的情況。
A=QR，正規正交矩陣Q不影響最小值，最小值取決於零餘部分R。
至於多解的情況，改為分解A的轉置矩陣Aᵀ=QR。
時間複雜度O(N³)。
但是計算量比NormalEquation少。
solveAx=b
min‖Ax-b‖²
min‖Qᵀ(Ax-b)‖²正規正交矩陣，變換後長度不變
min‖QᵀAx-Qᵀb‖²展開
min‖Rx-Qᵀb‖²QᵀA=QᵀQR=R
min‖[R₁x-Q₁ᵀb]‖²R=[R₁]區分出零，讓R₁是方陣
‖[0-Q₂ᵀb]‖[0]區分上段和下段
solveR₁x-Q₁ᵀb=0此式有唯一解，可為零（因此最小值是‖Q₂ᵀb‖²）
R₁x=Q₁ᵀb移項
x=R₁⁻¹Q₁ᵀb移項
數學公式（SingularValueDecompostion）
以下只證明無解的情況。
證明手法如出一轍。
至於多解的情況，改為分解A的轉置矩陣Aᵀ=UΣVᵀ。
時間複雜度O(N³)。
我不確定實務上是否比較快。
solveAx=b
min‖Ax-b‖²
min‖Uᵀ(Ax-b)‖²正規正交矩陣，變換後長度不變
min‖UᵀAx-Uᵀb‖²展開
min‖ΣVᵀx-Uᵀb‖²UᵀA=UᵀUΣVᵀ=ΣVᵀ
solveΣVᵀx-Uᵀb=0此式有唯一解，可為零
ΣVᵀx=Uᵀb移項
x=VΣ⁻¹Uᵀb移項
LinearLeastSquares:Optimization
TikhonovRegularization
一次方程組，無解、多解時，改為最佳化問題，以得到唯一解。
solveAx=b--->min‖Ax-b‖²[overdeterminedsystem]
solveAx=b--->min‖x‖²[underdeterminedsystem]
甚至利用Regularization，追加其他最佳化目標。
solveAx=b--->min‖Ax-b‖²+αf(x)(α≥0)
solveAx=b--->min‖x‖²+αf(x)(α≥0)
一次方程組，有許多式子和變數。
可能有其中一群變數與式子構成無解、另一群構成唯一解、剩下一群構成多解。
更有甚者，切割一些群結果無解變多解、整併某些群結果多解變無解。
無法釐清是無解、多解的時候，那就兩個一起上吧。
solveAx=b--->min‖Ax-b‖²+α‖x‖²(α≥0)
d/dx[‖Ax-b‖²+α‖x‖²]=0「一次微分等於零」的地方是極值、鞍點
二次函數、恆正，必得最小值
2AᵀAx-2Aᵀb+2αx=0展開
(AᵀA+αI)x=Aᵀb移項，左式即是AᵀA的對角線加上α
左式是實數對稱正定矩陣，有唯一解。
時間複雜度O(N³)。
HomogeneousLinearEquation
討論特例b=0的情況。
當b=0，則x=0，缺乏討論意義。
於是添加限制「x長度（的平方）為1」，增進討論意義。
solveAx=0--->min‖Ax‖²subjectto‖x‖²=1
min‖Ax‖²-λ(‖x‖²-1)Lagrangemultiplier
∂/∂x[‖Ax‖²-λ(‖x‖²-1)]=0「一次微分等於零」的地方是極值、鞍點
二次函數，必得極值
2AᵀAx-2λx=0展開
AᵀAx=λx移項，此即特徵向量的格式
答案是AᵀA的最小的特徵值的特徵向量！
AᵀA是對稱半正定矩陣，特徵值即是奇異值，特徵向量即是奇異向量。
AᵀA的特徵值即是A的奇異值，AᵀA的特徵向量即是A的右奇異向量。
LinearInequality
SystemofLinearInequalities
「一次不等式組」。
許多道一次不等式同時成立。
幾何意義是半空間交集，所有解形成一個凸多胞形，也可能形成開放區間、退化、空集合。
請見本站文件「Half-plane」。
目前沒有演算法。
大家更傾向討論「LinearProgramming」。
有人適度乘上負號，調整成Ax>b的格式，套用高斯消去法，但是不知道正不正確。
