Sparse Linear Equation
Sparse Linear Equation
「稀疏一次方程式」。一次方程式，給定矩陣是稀疏矩陣。
solve Ax = b where ‖A‖₀ ≤ k
專著《Direct Methods for Sparse Linear Systems》。
Matrix Splitting
「分裂」。Ax = b變成(M + N)x = b。原本矩陣分裂成兩個矩陣相加，使得M和N清爽暢快，容易計算。
(M + N)x = b M⁻¹(M + N)x = M⁻¹b x + M⁻¹Nx = M⁻¹b x = -M⁻¹Nx + M⁻¹b
M必須是可逆矩陣（擁有反函數）。
得到遞推公式xnext = -(M⁻¹N)x + (M⁻¹b)。不保證收斂。
Matrix Preconditioning
「預條件」。分為左右兩種。
左：Ax = b變成MAx = Mb。等號兩邊同乘矩陣M，使得MA清爽暢快，容易計算。先算MA與Mb，再解(MA)x = (Mb)。
右：Ax = b變成AMz = b。向量x代換成Mz，使得AM清爽暢快，容易計算。先解(AM)z = b，再求x = Mz。
M必須是可逆矩陣（擁有反函數），讓解的數量不變。
左右可以個別使用，也可以一齊使用。
Sparse Linear Equation: Permutation
Reverse Cuthill–McKee Algorithm
細節請見維基百科。此處僅簡介。
Band Matrix：元素盡量集中於對角線。
解方程組，可以自由交換橫條與直條。
交換橫條與直條，調整成帶狀矩陣。減少cache miss。
Bandwidth：每個橫條，非零區間寬度，最大值。
帶寬最窄的矩陣，很不幸的，NP-complete。大家只好採用近似演算法，快速得到結果，例如Reverse Cuthill–McKee Algorithm。
Dulmage–Mendelsohn Decomposition
細節請見課程講義。此處僅簡介。
Block Triangular Matrix：分塊視作元件，形成三角矩陣。
交換橫條與直條，調整成分塊三角矩陣。減少求解時間。
Structural Rank：最大二分匹配的邊數。
矩陣視作二分圖，找到二分圖最大匹配，重新整理成三個部分：overdetermined system、underdetermined system、consistent system。
Birkhoff-von Neumann Decomposition
細節請見課程講義。此處僅簡介。
Doubly Stochastic Matrix：各橫條總和皆為1、各直條總和皆為1。橫條暨直條，故取名Doubly，總和為一宛如機率，故取名Stochastic。
雙隨機矩陣可以分解成排列矩陣們的加權平均數。權重總和為1，權重都是0到1。通常有許多種。
數學式子形成Splitting的形式。
A = w₁ P₁ + w₂ P₂ + ... + wₖ Pₖ
排列矩陣數量最少的BvN分解，很不幸的，NP-complete。大家只好採用近似演算法，快速得到結果。
排列矩陣數量最多需要(N-1)²+1個。藉此得到時間複雜度。
Sparse Linear Equation: Relaxation
Jacobi Iteration
細節請見本站文件「Linear Equation」。此處僅簡介。
鬆弛法的本質是Splitting。Ax = b變成(D+L+U)x = b。遞推公式xnext = -D⁻¹(L+U)x + D⁻¹b。
鬆弛法僅適用於對角優勢矩陣。優點是容易實作，缺點是末期收斂速度慢。
Sparse Linear Equation: Decomposition
LU Decomposition
Cholesky Decomposition
細節請見課程講義。此處僅簡介。
高斯消去法，難以保持稀疏。LU分解，有辦法盡量保持稀疏。
一般矩陣：LU分解。對稱正定矩陣：Cholesky分解。
A = LU LU Decomposition A = LLᵀ Cholesky Decomposition (unique solution)
解方程組，可以自由交換橫條與直條。藉由交換，改變矩陣，改變LU分解。我們希望找到一種交換方式，讓L和U盡量稀疏。
寫成數學式子，就是套用排列矩陣。PA交換橫條，APᵀ交換直條。注意到，交換直條APᵀ，也需要交換變數Pb。PAPᵀ同時交換橫條直條，以保持對稱。最後形成Preconditioner的形式。
PAQ = Qb preconditioning PAPᵀ = Pb preconditioning of symmetric matrix
LU分解較複雜，Cholesky分解較簡單，此處僅討論後者。
Cholesky分解可以用圖論來解釋。已知原圖A，求得新圖L+Lᵀ。每回合消除一點i。消除之前若存在邊ix、邊iy，消除之後將出現邊xy。
已知原圖A，找到一種消除順序，盡量不增加邊，讓新圖L+Lᵀ盡量稀疏。最佳消除順序稱作Minimum Elimination Ordering，最少需要增加的邊稱作Minimum Fill-in。
最佳消除順序，很不幸的，NP-complete。大家只好採用近似演算法，快速得到結果，例如Minimum Degree Algorithm。
順帶一提，只增加0條邊，消除順序正是Perfect Elimination Ordering，矩陣正是Chordal Graph。
消除順序，畫成一棵樹，稱作Elimination Tree。它是新圖的Transitive Reduction，也是新圖的DFS Tree。在圖論當中，一旦有了樹，就能變出許多把戲。
Cholesky分解的計算過程，就是Elimination Tree的拓樸順序，從葉往根。平行演算法：Multifrontal、Supernodal。
http://buttari.perso.enseeiht.fr/docs/ALC_Toulouse.pdf http://faculty.cse.tamu.edu/davis/research_files/siam06_short.pdf
知名函式庫如CHOLMOD、Intel MKL、PARDISO。
Incomplete LU Decomposition
Incomplete Cholesky Decomposition
故意不增加邊，維持0。計算結果錯誤，但是算得快。
Sparse Linear Equation: Optimization
Preconditioned Biconjugate Gradient Method
Preconditioned Conjugate Gradient Method
儘管Incomplete LU Decomposition結果不正確，但是可以拿來當作Preconditioner。
Ax = b變成(LU)⁻¹Ax = (LU)⁻¹b。(LU)⁻¹A差不多是I，導致收斂步數相當少。
順帶一提，反矩陣近似於原矩陣、又很稀疏的Preconditioner，有些文獻稱作Sparse Approximate Inverse Preconditioner。
一般矩陣：ILU配PBiCG。對稱正定矩陣：IChol配PCG。
IChol和PCG巧妙避開O(N³)。實務上最快的演算法。
Multigrid Preconditioned Conjugate Gradient Method
追加兩個外掛，硬是催出效能：Precondition調整矩陣數值、Multigrid調整矩陣大小。
Precondition是以等量乘法公理，讓原本矩陣乘上自訂矩陣，讓原本矩陣變漂亮、容易求解。
Multigrid就是Scaling Method的意思。原本矩陣縮放成各種大小（各種解析度），從小到大，依序計算，整合結果。
詳細內容我沒有學會。大家自己看著辦吧。
可以直接使用MATLAB的pcg()。
Sparse Linear Least Squares
Sparse Linear Least Squares
無解、多解時，改為找到平方誤差最小的解。
min ‖Ax - b‖₂² where ‖A‖₀ ≤ k [no solution] min ‖x‖₂² where ‖A‖₀ ≤ k [infinitely many solutions]
專著《Iterative Methods for Sparse Linear Systems》。
Matrix Preconditioning: Aᵀ
Aᵀ的功能是投影降維。無解、多解時，利用左右預條件矩陣Aᵀ調整維度，以得到唯一解。
無解：左預條件。輸出降維，左邊乘以Aᵀ。
多解：右預條件。輸入升維，右邊除以Aᵀ。
| equation | optimization
------------| ------------------------| -------------------------
0 solution | solve AᵀAx = Aᵀb (over) | min ‖Ax - b‖² (residual)
∞ solutions | solve AAᵀz = b (under)| min ‖x - x*‖² (error)
------------| ------------------------| -------------------------
1 solution | solve Ax = b (well) | ‖Ax - b‖² = ‖x - x*‖² = 0
其中x*是最佳解位置。多解時，x* = 0。【尚待確認】
線性函數兩大問題：求解、求特徵向量暨特徵值
一般函數，求根、求不動點、求特徵點、求解，四者等價。
線性函數，四者難度不同。大家討論求解（涵蓋求根）、求特徵點（涵蓋不動點）。兩者各有演算法，求解較易，求特徵點較難。
find solution Ax = b find eigenvector/value Ax = λx
順帶一提，特徵向量暨特徵值可用來求解。
線性函數兩大進階問題：求最小平方解、求奇異向量暨奇異值
左右預條件矩陣Aᵀ。
[left preconditioner Aᵀ] overdetermined least squares Aᵀ(Ax) = Aᵀ(b) right singular vector/value Aᵀ(Ax) = Aᵀ(σx) [right preconditioner Aᵀ] underdetermined least squares A(Aᵀz) = b where x = Aᵀz left singular vector/value A(Aᵀz) = σ(Aᵀz) where x = Aᵀz
順帶一提，奇異向量暨奇異值可用來求最小平方解。
Sparse Linear Least Squares: Decomposition
演算法（CGNR）
Normal Equation + Conjugate Gradient Method。
CGNR: left preconditioner AᵀAx = Aᵀb CGNE: right preconditioner AAᵀz = b
演算法（Generalized Minimal Residual Method）（GMRES）
只能用於方陣。
一、變成幾乎上三角矩陣。二、變成上三角矩陣。三、求解。
一、一般矩陣的幾乎上三角分解：Arnoldi Iteration，初始向量設為b。原矩陣A簡化成幾乎上三角矩陣H，原向量b簡化成只有第一個元素有值的向量‖b‖e₁。
solve Ax = b => min ‖Ax - b‖² => min ‖QHQᵀx - b‖² Hessenberg decomposition => min ‖HQᵀx - Qᵀb‖² multiply Qᵀ and norm remain unchanged => min ‖Hz - Qᵀb‖² let Qᵀx = z => solve Hz = Qᵀb
Qᵀb = [q₁∙b₁ q₂∙b₂ ... qₙ∙bₙ]ᵀ
= [‖b‖ 0 0 0 0 ...]ᵀ let x₀ = b and q₁ = x₀/‖x₀‖
= ‖b‖e₁ let e₁ = [1 0 0 0 0 ...]ᵀ
二、幾乎上三角矩陣的QR分解：適合使用Givens Rotation。
solve Ax = b A := H and b := ‖b‖e₁ => min ‖Ax - b‖² => min ‖QRx - b‖² QR decomposition => min ‖Rx - Qᵀb‖² multiply Qᵀ and norm remain unchanged => solve Rx = Qᵀb
三、上三角矩陣求解：back substitution。
high school math
順帶一提，步驟一也可以用Galerkin Method來解釋。
solve Ax = b => solve Qᵀ(Ax - b) = 0 Galerkin method => solve QᵀAx = Qᵀb => solve QᵀAQz = Qᵀb let x = Qz => solve Hz = Qᵀb
順帶一提，步驟一也可以將初始向量設為r₀ = b - Ax₀，其中x₀是任意值。原向量b簡化成只有第一個元素有值的向量‖r₀‖e₁。
solve Ax = b => solve Qᵀ(Ax - b) = 0 Galerkin method => solve Qᵀ(Ax₀ + AQz - b) = 0 let x = x₀ + Qz => solve QᵀAQz + Qᵀ(Ax₀ - b) = 0 => solve QᵀAQz - Qᵀr₀ = 0 let r₀ = b - Ax₀ => solve QᵀAQz = Qᵀr₀ => solve Hz = Qᵀr₀
演算法（Minimal Residual Method）（MINRES）
只能用於對稱方陣。Arnoldi Iteration改成Lanczos Iteration。
Lanczos Iteration的時間複雜度分為兩部分：矩陣求值總共O(N³)、向量扳正總共O(N²)。稀疏矩陣的情況下，矩陣求值需時變少，得以發揮長處。
演算法（LSMR）
Normal Equation + Minimal Residual Method。
solve AᵀAx = Aᵀb normal equation (AᵀA is symmetric) => solve QHQᵀx = Aᵀb Hessenberg decomposition => solve HQᵀx = QᵀAᵀb multiply Qᵀ and norm remain unchanged => solve Hz = QᵀAᵀb let Qᵀx = z => solve Hz = (AQ)ᵀb
似乎是實務上最快的演算法。【尚待確認】
演算法（LSQR）
一、變成雙對角線矩陣。二、變成上三角矩陣。三、求解。
solve Ax = b
=> min ‖Ax - b‖²
=> min ‖VBUᵀx - b‖² Golub–Kahan bidiagonalization
=> min ‖BUᵀx - Vᵀb‖² multiply Vᵀ and norm remain unchanged
=> min ‖Bz - Vᵀb‖² let Uᵀx = z
=> min ‖ ⎡ B₁z - Vᵀb ⎤ ‖² let B = ⎡ B₁ ⎤ , B₁ is square matrix
‖ ⎣ 0 - Vᵀb ⎦ ‖ ⎣ O ⎦
=> min ‖B₁z - Vᵀb‖²
=> solve B₁z = Vᵀb