Sparse Matrix

Sparse Matrix

「稀疏矩陣」。矩陣元素幾乎都是零。非零元素相當稀少。

矩陣的資料結構=圖的資料結構。稠密矩陣適合Adjacency Matrix,稀疏矩陣適合Adjacency Lists。

實務上採用Adjacency Lists與Edge List的融合版本,化作三條陣列。MATLAB即採用此方式。

Sparse Eigenvalue

Sparse Eigenvalue

經典的矩陣運算,諸如反矩陣、QR分解、特徵分解,時間複雜度O(N³)以上,速度太慢不實用。因此數學家旁敲側擊,鑽研稀疏矩陣,精簡計算步驟,以投入實用。

專著《Numerical Methods for Large Eigenvalue Problems》

幾乎上三角分解演算法(Restarted Arnoldi Iteration)

一口氣求得大量特徵值,其中一種策略是:先變成幾乎上三角矩陣、再變成上三角矩陣。

幾乎上三角矩陣分解,三種演算法Arnoldi Iteration、Householder Reflection、Jacobi Rotation,其中只有Arnoldi Iteration適用於稀疏矩陣。

Arnoldi Iteration的時間複雜度,分成兩個部分。

一、矩陣求值:稀疏矩陣的情況下,需時較少。

二、向量扳正:無法保持稀疏,需時相同。

數學家想到了歪招:restart。被扳正的向量,每回合越來越多,計算量越來越大。因此每隔幾回合就砍掉重練,初始向量換成新的再來一遍,直到總共湊滿N個向量。

restart也能解決「初始向量不足N種特徵分量,無法生成整個空間」的情況。

Large Scale Matrix

Large Scale Matrix

「大型矩陣」。矩陣尺寸非常非常大。必須煩惱如何存取。

http://tmtacm.blogspot.com/2016/01/blog-post_26.html

Large Scale Eigenvalue

Large Scale Eigenvalue

採用平行計算。

特徵向量暨特徵值修正演算法(Olsen's Method)

特徵向量暨特徵值,寫成方程組。

{ Ax = λx
{ ‖x‖² = 1

移項整理,求解改成求根。

{ (A-λI)x = 0
{ ½xᵀx - 1 = 0     為了方便微分,特徵向量長度改成2。

牛頓法求根。

xnext = x - F′(x)ᵀ⁺ F(x)

[ xnext ] = [ x ] - [ df₁/dx df₂/dx ]ᵀ⁻¹[ f₁(x,λ) ]
[ λnext ]   [ λ ]   [ df₁/dλ df₂/dλ ]   [ f₂(x,λ) ]

[ xnext ] = [ x ] - [ (A-λI)ᵀ x ]ᵀ⁻¹[ r ]
[ λnext ]   [ λ ]   [   -xᵀ   0 ]   [ 0 ]

[ xnext ] = [ x ] - [ A-λI -x ]⁻¹[ r ]   where r = (A-λI)x = Ax - λx
[ λnext ]   [ λ ]   [  xᵀ   0 ]  [ 0 ]

簡易方式,結果相同。

A(x + Δx) = (λ + Δλ)(x + Δx)

(A-λI)Δx - xΔλ = r     where r = Ax - λx and ΔλΔx ≈ 0

[ A-λI -x ] [ Δx ] = [ r ]     where Δx = xnext - x
[  xᵀ   0 ] [ Δλ ]   [ 0 ]           Δλ = λnext - λ

每回合都要矩陣求解,時間複雜度相當高,似乎沒有用處。

特徵向量暨特徵值修正演算法(Rayleigh Quotient Iteration)

Rayleigh Quotient Iteration是牛頓法的變種。

A(x + Δx) = (λ + Δλ)(x + Δx)
A(x + Δx) = λ(x + Δx) + Δλ(x + Δx)
(A - λI)(x + Δx) = Δλ(x + Δx)
(A - λI)(x + Δx) = Δλx                 where ΔλΔx ≈ 0
(x + Δx) = (A - λI)⁻¹Δλx
xnext = Δλ(A - λI)⁻¹x       let λ = xᵀAx / xᵀx is Rayleigh quotient

「向量長度縮放倍率Δλ」改成了「向量長度縮放為1」。

「λnext = λ + Δλ」改成了「λnext是虛擬特徵值」。

每回合都要矩陣求解,時間複雜度相當高,似乎沒有用處。

特徵向量暨特徵值修正演算法(Jacobi-Davidson Method)

追加垂直限制。

A(x + Δx) = (λ + Δλ)(x + Δx)    where x ⟂ Δx

移項整理。

(A-λI)Δx - Δλx = r     where x ⟂ Δx, r = Ax - λx

利用投影矩陣,消掉-Δλx。

let P = I - xxᵀ     such that Px = 0 and Pr = r
P((A-λI)Δx - Δλx) = Pr
P(A-λI)Δx - PΔλx = Pr
P(A-λI)Δx = r

補上右邊,成為相似變換,保留特徵值。

Pᵀ(A-λI)PΔx = r     (P = Pᵀ and PΔx = Δx)

問題變成一次方程組求解。求得特徵向量誤差Δx。

BΔx = r     where B = Pᵀ(A-λI)P and r = (A-λI)x = Ax - λx

每回合都要矩陣求解,時間複雜度相當高,似乎沒有用處。