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)
解方程組，可以自由交換橫條與直條。藉由交換，改變矩陣A，進而改變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的意思。原本矩陣縮放成各種大小（各種解析度），從小到大（從粗到細），依序計算，整合結果。
為了縮短計算時間，分成三種計算順序：V-cycle、F-cycle、W-cycle。前者省時但答案不精確，後者耗時但答案較精確。
示意圖請見谷歌搜尋：
可以直接使用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
演算法（conjugate gradient method on the normal equation）（CGNE）
normal equation + conjugate gradient method。
overdetermined least squares (left preconditioner): AᵀAx = Aᵀb underdetermined least squares (right preconditioner): AAᵀz = b
改良版本CGLS將整件事情視作preconditioned conjugate gradient method。
演算法（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 Az = b A := H and b := ‖b‖e₁ => min ‖Az - b‖² => min ‖QRz - b‖² QR decomposition => min ‖Rz - Qᵀb‖² multiply Qᵀ and norm remain unchanged => solve Rz = Qᵀb
三、上三角矩陣求解：back substitution。
solve Rz = y y := Qᵀb => z = high school math
四、正規正交矩陣求解：matrix–vector multiplication。
solve Qᵀx = z => x = Qz since Q⁻¹ = Qᵀ
順帶一提，步驟一也可以用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²)。稀疏矩陣的情況下，矩陣求值需時變少，得以發揮長處。
演算法（MINRES on the normal equation）
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₁ ⎤ and B₁ is square
‖⎣ 0 - Vᵀb ⎦‖ ⎣ O ⎦
=> min ‖B₁z - Vᵀb‖²
=> solve B₁z = Vᵀb
改良版本LSMR微調了QR分解的計算步驟。