sparse matrix

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

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

實務上採用adjacency lists與edge list的融合版本，化作三條陣列，稱作compressed sparse row。MATLAB即採用此方式。

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 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

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