sparse matrix

sparse matrix

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

矩陣的資料結構=圖的資料結構。稠密矩陣適合adjacency matrix,稀疏矩陣適合adjacency lists。

實務上採用adjacency lists與edge list的融合版本,化作三條陣列,稱作compressed sparse row。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

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

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

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