Sparse Linear Equation(Under Construction!)

Sparse Linear Equation

「稀疏一次方程式」。一次方程式,給定矩陣是稀疏矩陣。

solve Ax = b   where ‖A‖₀ ≤ k

以下將介紹基礎知識、相關問題。

Matrix Splitting

「矩陣分裂」。A = A₁ + A₂ + ...,原本矩陣分裂成多個矩陣相加,使得A₁ A₂ ...清爽暢快,容易計算。

移項法則,換個角度來思考,即是A = D + L + U。

Preconditioner

「預條件」。Ax = b變成MAx = Mb,等號兩邊同乘矩陣M,使得MA清爽暢快,容易計算。

M必須是可逆矩陣(擁有反函數),讓解的數量不變。

移項法則,從某種程度上而言,即是等號兩邊同乘矩陣D⁻¹。

LU Decomposition(Under Construction!)

Sparse Approximate Inverse Preconditioner

http://www.mathcs.emory.edu/~benzi/Web_papers/comp.pdf

請見專著《Direct Methods for Sparse Linear Systems》

Sparse LU Decomposition

Sparse Cholesky Decomposition

http://www.cs.tau.ac.il/~stoledo/taucs/
Supernodal Method
Multifrontal Method
CHOLMOD

Incomplete LU Decomposition

Incomplete Cholesky Decomposition

https://rosettacode.org/wiki/Cholesky_decomposition
https://www.geeksforgeeks.org/cholesky-decomposition-matrix-decomposition/

Relaxation(Under Construction!)

QR Decomposition(Under Construction!)

Sparse Linear Least Squares

請見專著《Iterative Methods for Sparse Linear Systems》

Krylov Matrix

唸做kree-luv而非cry-love。

原矩陣必須N維,Krylov Matrix才是N維。

     [  |   |        |   ]
K  = [ A⁰b A¹b ... Aⁿ⁻¹b ] = [y₀ y₁ ... yₙ₋₁]
     [  |   |        |   ]

     [  |   |       |  ]
AK = [ A¹b A²b ... Aⁿb ]   = [y₁ y₂ ... yₙ]
     [  |   |       |  ]

            [  |  |      |  | ]
AK = KC = K [ e₁ e₂ ... eₙ -c ]     where c = -K⁻¹Aⁿy₀
            [  |  |      |  | ]

            [ 0 0 ... 0 -c₀   ]
            [ 1 0 ... 0 -c₁   ]
K⁻¹AK = C = [ 0 1 ... 0 -c₂   ]
            [ : :     :  :    ]
            [ 0 0 ... 1 -cₙ₋₁ ]
enlarged Krylov subspace
https://who.rocq.inria.fr/Laura.Grigori/
https://people.eecs.berkeley.edu/~demmel/
rational Krylov subspace
http://guettel.com/rktoolbox/examples/html/index.html

Similarity Transformation

Q⁻¹AQ = B。原本矩陣A,相似變換矩陣Q,相似矩陣B。

相似矩陣特徵值相同。相似矩陣可用來求特徵值。

Rayleigh-Ritz Method:任意正規正交矩陣拿來當作Q。
        QR Iteration:A的QR分解拿來當作Q。計算公式可以精簡成B=RQ。
   Arnoldi Iteration:K的QR分解拿來當作Q。B會是幾乎上三角矩陣。

演算法(Arnoldi Iteration)

http://www.cs.cornell.edu/courses/cs6220/2017fa/CS6220_Lecture8.pdf
                             [ R₁₁ R₁₂ R₁₃ .. .. ]
            [  |  |      | ] [  0  R₂₂ R₂₃ .. .. ]
Kₙ = QₙRₙ = [ q₁ q₂ ... qₙ ] [  0   0  R₃₃ .. .. ]
            [  |  |      | ] [  :   :   :  \    ]
                             [  :   :   :     Rₙₙ]

                                             [ H₁₁ H₁₂ H₁₃ .. .. ]
      [  |   |       | ]   [  |  |      |  ] [ H₂₁ H₂₂ H₂₃ .. .. ]
AQₙ = [ Aq₁ Aq₂ ... Aqₙ ] = [ q₁ q₂ ... qₙ₊₁] [  0  H₃₂ H₃₃ .. .. ] = Qₙ₊₁Hₙ
      [  |   |       | ]   [  |  |      |  ] [  0   0  H₃₄ .. .. ]
                                             [  :   :   :  \    ]
                                             [  :   :   :    Hₙ₊₁ₙ]
Aqₙ = H₁ₙq₁ + H₂ₙq₂ + ... + Hₙ₊₁ₙqₙ₊₁
qᵢᵀAqₙ = Hᵢₙ
q₁ = u/‖u‖
for n = 1,2,...
   v := Hₙ₊₁ₙqₙ₊₁ = Aqₙ - H₁ₙq₁ - H₂ₙq₂ - ... - Hₙₙqₙ
   Hₙ₊₁ := ‖v‖
   qₙ₊₁ := v/‖v‖

K = QR。

Krylov矩陣,做QR分解,演算法是Gram-Schmdit正交化。

H = Q⁻¹AQ。

不算R,改算H。H是幾乎上三角矩陣。

H = Q⁻¹AQ = RCR⁻¹。

H = QᵀAQ是二次型。

A Qn = Qn+1 Hn。

演算法(Lanczos Iteration)

A是對稱矩陣,計算公式可精簡。

對稱矩陣改用Wilkinson Shift。

Aqₙ = Hₙ₋₁qₙ + Hₙₙqₙ + Hₙ₊₁ₙqₙ₊₁

Linear Least Squares

QR Decomposition ---> QR Iteration。

AK:GMRES/MINRES,K:FOM/SYMMLQ。

演算法(Generalized Minimum Residual Method)(GMRES)

http://www.math.ubbcluj.ro/~tradu/nlalgslides/Krylovsmeth.pdf

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

一邊Arnoldi Iteration,一邊求解。以Givens Rotation加速求解。

  min ‖Ax - b‖²
= min ‖A (Kn c) - b‖²       let x = Kn c
= min ‖A (Qn z) - b‖²       let x = Qn z   (Kn = Qn Rn)
= min ‖Qn+1 Hn z - b‖²
= min ‖Hn z - Qn+1⁻¹ b‖²    multiplication of orthonormal matrix
= min ‖Hn z - Qn+1ᵀ b‖²     Qn+1⁻¹ will not change norm
= min ‖Hn z - ‖b‖ e₁‖²
= min ‖Hn z - [‖b‖ 0 0 0 ...]ᵀ‖²

額外介紹兩個小技巧:

shift:暫時調整特徵值,加快收斂速度。

restart:Q越長越大,記憶體不足,於是每隔幾回合就砍掉重練,Krylov Matrix種子換成新的。

演算法(Minimum Residual Method)(MINRES)

對稱矩陣。

Optimization(Under Construction!)

Linear Approximation【尚無正式名稱】

pseudoinverse = projection = least squares = linear regression
= dot product of arbitrary basis

linear approximation
= dot product of othronormal basis
galerkin method: orthonormal basis dot
conjugated gradient method: A-orthonormal residual dot
linear feedback system: time series dot (fixedpoint)
companion matrix: x x^2 x^3 dot (differential equation eigenfunction)
taylor expansion: order-1 to order-n gradient dot
kalman filter: residual dot
https://www.encyclopediaofmath.org/index.php/Galerkin_method
正規正交基底的線性組合的線性變換,分別對各個基底點積。
Ax = A(c1q1 + c2q2 + ...)
q*Ax = c q*Aq   when dot(qi,qj) = 0

點積結果是二次型。正規化之後得到虛擬特徵值,乘上c倍。
c q*Aq/q*q

https://en.wikipedia.org/wiki/QR_decomposition
QR分解之GS正交化 Q當作正規正交基底 R是點積
dot(q,a) = q*a
[galerkin method solve linear equation]
xk = x0 + subspace = x0 + (α0 b + α1 Ab + α2 A²b)
such that dot(b - Axk, q) = 0   餘數通通通過基底q的測試
如果是 linear least squares  應該要改成 min dot(b - Axk, q)
xₜ₊₁ = xₜ - dot(J, F) / dot(J, J)
視作投影公式:F(x)投影到J。
向量投影到梯度,再求得向量到梯度的差距。
can we use them? when use the linearization
1. fixpoint iteration (linear equation)
2. power iteration (linear eigenvalue)
n->n求根   其實就是求特徵值? power iteration?
           companion matrix? Cayley-Hamilton Theorem?

演算法(Biconjugate Gradient Method)

正定矩陣。

演算法(Conjugate Gradient Method)

對稱正定矩陣。

http://www.cse.iitd.ac.in/~narain/courses/col865i1819/slides/06-equation-solving-optimization.pdf

solve Ax = b
span {b}         : assume x = α0 b                 , find best α0
span {b, Ab}     : assume x = α0 b + α1 Ab         , find best α0, α1
span {b, Ab, A²b}: assume x = α0 b + α1 Ab + α2 A²b, find best α0, α1, α2
......

Sparse Linear Least Squares(Under Construction!)

Sparse Linear Least Squares

無解、多解時,改為找到平方誤差最小的解。

min ‖Ax - b‖₂²                  where ‖A‖₀ ≤ k   [overdetermined system]
min ‖x‖₂²   subject to Ax = b   where ‖A‖₀ ≤ k   [underdetermined system]

演算法(LSQR)

https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.sparse.linalg.lsqr.html#scipy.sparse.linalg.lsqr

演算法(LSMR)

https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.sparse.linalg.lsmr.html#scipy.sparse.linalg.lsmr