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

https://en.wikipedia.org/wiki/Band_matrix

細節請見維基百科。此處僅簡介。

Band Matrix:元素盡量集中於對角線。

解方程組,可以自由交換橫條與直條。

交換橫條與直條,調整成帶狀矩陣。減少cache miss。

Bandwidth:每個橫條,非零區間寬度,最大值。

帶寬最窄的矩陣,很不幸的,NP-complete。大家只好採用近似演算法,快速得到結果,例如Reverse Cuthill-McKee Algorithm。

Dulmage-Mendelsohn Decomposition

http://perso.ens-lyon.fr/bora.ucar/CR02/

細節請見課程講義。此處僅簡介。

Block Triangular Matrix:區塊視作元件,形成三角矩陣。

交換橫條與直條,調整成區塊三角矩陣。減少求解時間。

Structural Rank:最大二分匹配的邊數。

矩陣視作二分圖,找到二分圖最大匹配,重新整理成三個部分:overdetermined system、underdetermined system、consistent system。

Birkhoff-von Neumann Decomposition

http://perso.ens-lyon.fr/bora.ucar/CR02/

細節請見課程講義。此處僅簡介。

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

https://en.wikipedia.org/wiki/Matrix_splitting#Matrix_iterative_methods

細節請見本站文件「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

http://perso.ens-lyon.fr/bora.ucar/CR02/

細節請見課程講義。此處僅簡介。

高斯消去法,難以保持稀疏。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 

知名函式庫如CHOLMODIntel MKLPARDISO

Incomplete LU Decomposition
Incomplete Cholesky Decomposition

故意不增加邊,維持0。計算結果錯誤,但是算得快。

Sparse Linear Equation: Optimization

Preconditioned Biconjugate Gradient Method
Preconditioned Conjugate Gradient Method

https://scicomp.stackexchange.com/questions/26862/

儘管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³)。實務上最快的演算法。

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 solution   Aᵀ(Ax) = Aᵀ(b)
right singular vector/value             Aᵀ(Ax) = Aᵀ(σx)

[right preconditioner Aᵀ]
underdetermined least squares solution  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 Minimum 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 v₁ = b and q₁ = v₁/‖v₁‖
    = ‖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₀

演算法(Minimum Residual Method)(MINRES)

只能用於對稱方陣。Arnoldi Iteration改成Lanczos Iteration。

Lanczos Iteration的時間複雜度分為兩部分:矩陣求值總共O(N³)、向量扳正總共O(N²)。稀疏矩陣的情況下,矩陣求值需時變少,得以發揮長處。

演算法(LSMR)

Normal Equation + Minimum 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