數值分析系列二：進階特徵值方法與應用

特徵值問題是數值分析的核心，應用由結構工程延伸至量子力學。簡潔的方程

Ax=\lambda x

背後包含穩定性、計算成本、矩陣結構與尺度等問題。本文由基本迭代走到大型稀疏方法，再概覽工程、量子化學、機器學習及網絡分析應用。

基本特徵值方法

冪法

冪法反覆計算

x_{k+1}=\frac{Ax_k}{\|Ax_k\|}, \qquad \lambda_k=\frac{x_k^TAx_k}{x_k^Tx_k}.

若矩陣可對角化、主導特徵值在模長上唯一，而且初始向量在其特徵向量方向有非零分量，向量方向通常以約

\left|\frac{\lambda_2}{\lambda_1}\right|^k

的比例收斂。Rayleigh quotient 的特徵值誤差對對稱矩陣往往可以有更快的局部階，因此不能把同一界線不加條件地套用到所有 $\lambda_k$ 。

位移反迭代

要尋找最接近指定 $\mu$ 的特徵值，可求解

(A-\mu I)y_{k+1}=x_k, \qquad x_{k+1}=\frac{y_{k+1}}{\|y_{k+1}\|}.

因為每一步使用同一個位移矩陣，可以先作一次 LU 分解：

def inverse_iteration(A, mu, max_iter=100, tol=1e-10):
    n = A.shape[0]
    x = np.random.rand(n)
    x /= np.linalg.norm(x)
    I = np.eye(n)
    lu, piv = scipy.linalg.lu_factor(A - mu*I)

    for k in range(max_iter):
        y = scipy.linalg.lu_solve((lu, piv), x)
        x_new = y / np.linalg.norm(y)

        # Eigenvectors are unchanged by a sign flip.
        residual = min(
            np.linalg.norm(x_new - x),
            np.linalg.norm(x_new + x)
        )
        x = x_new
        if residual < tol:
            eigenvalue = (x.T @ A @ x) / (x.T @ x)
            return eigenvalue, x

    return (x.T @ A @ x) / (x.T @ x), x

原始示例以 mu + 1 / ||y|| 估計特徵值，但向量符號與正規化會令這個式子不穩妥；直接使用 Rayleigh quotient 較清楚。若 $\mu$ 極接近特徵值，線性系統本身會病態，求解精度亦需檢查。

QR 演算法

實際密集矩陣演算法通常先把 $A$ 約化成上 Hessenberg 形式：

A=QHQ^T.

其後作帶位移 QR 迭代：

H_k-\mu_kI=Q_kR_k, \qquad H_{k+1}=R_kQ_k+\mu_kI.

這是一個相似變換，所以保持特徵值。當下次對角元素足夠小，矩陣逐步分裂成較小區塊；實數矩陣最後形成 real Schur form，由 $1\times1$ 與 $2\times2$ 區塊讀出實數及共軛複數特徵值。

QR 的全局與局部收斂率取決於矩陣類型、特徵值分離、位移策略及 deflation；不能簡化成只由 $H_0$ 的譜半徑控制的一條普遍幾何界線。

兩個定位與擾動定理

Gershgorin 圓盤定理

對 $A=[a_{ij}]\in\mathbb C^{n\times n}$ ，每個特徵值至少位於一個圓盤：

|z-a_{ii}|\leq\sum_{j\ne i}|a_{ij}|.

它提供快速而便宜的譜位置界線，也能協助選擇位移與檢查不合理輸出。

Bauer–Fike 定理

若 $A=X\Lambda X^{-1}$ 可對角化，對任意擾動 $E$ ， $A+E$ 的每個特徵值 $z$ 都滿足

\min_{\lambda\in\sigma(A)}|z-\lambda| \leq\kappa(X)\|E\|.

$\kappa(X)$ 是特徵向量矩陣的條件數。正常矩陣可選酉矩陣 $X$ ，條件數為 1；高度非正常矩陣則可能令很小的 $E$ 造成很大譜移動。這個定理談的是 $A+E$ 的特徵值，而不是把近似特徵值直接放入一個自造的對角殘差矩陣。

Deflation 與位移

實作常在

|h_{i+1,i}| \leq\epsilon\left(|h_{ii}|+|h_{i+1,i+1}|\right)

時把次對角項視為可忽略，再 deflate 已收斂區塊。常見位移包括 Rayleigh quotient、Wilkinson shift 與處理停滯的 exceptional shifts。

對實對稱矩陣尾端 $2\times2$ 區塊，Wilkinson shift 可寫成：

def wilkinson_shift(H):
    a = H[-2, -2]
    b = H[-1, -1]
    c = H[-2, -1]
    d = (a - b) / 2.0

    if d == 0:
        return b - abs(c)
    sign = 1 if d > 0 else -1
    return b - c*c / (d + sign*np.sqrt(d*d + c*c))

這個簡式假設尾端對稱；一般非對稱 Hessenberg QR 需要處理 $2\times2$ 區塊和雙位移。

利用矩陣結構

對稱矩陣與 Lanczos

實對稱矩陣有實特徵值及正交特徵向量。Lanczos 以三項遞迴建立 Krylov 子空間，得到小型三對角投影：

def lanczos_iteration(A, v0, k):
    n = A.shape[0]
    V = np.zeros((n, k + 1))
    T = np.zeros((k + 1, k))
    V[:, 0] = v0 / np.linalg.norm(v0)

    for j in range(k):
        w = A @ V[:, j]
        if j > 0:
            w -= T[j - 1, j - 1] * V[:, j - 1]

        T[j, j] = np.dot(w, V[:, j])
        w -= T[j, j] * V[:, j]

        # Full reorthogonalization for this teaching version.
        for i in range(j):
            w -= np.dot(w, V[:, i]) * V[:, i]

        beta = np.linalg.norm(w)
        if j + 1 < k:
            T[j + 1, j] = beta
        if beta < 1e-14:
            return V[:, :j + 1], T[:j + 1, :j + 1]
        V[:, j + 1] = w / beta

    return V[:, :k], T[:k, :k]

原始短碼的 T 索引把前一步係數放在錯誤位置。以上仍只是教學版本；成熟實作要處理重正交化、鎖定已收斂 Ritz vectors 與稀疏 matvec。

非對稱稀疏矩陣與 Arnoldi

Arnoldi 建立正交基底 $V$ 及上 Hessenberg 投影 $H$ ：

def arnoldi_iteration(A, v0, k):
    n = A.shape[0]
    V = np.zeros((n, k + 1))
    H = np.zeros((k + 1, k))
    V[:, 0] = v0 / np.linalg.norm(v0)

    for j in range(k):
        w = A @ V[:, j]
        for i in range(j + 1):
            H[i, j] = np.dot(w, V[:, i])
            w -= H[i, j] * V[:, i]

        H[j + 1, j] = np.linalg.norm(w)
        if H[j + 1, j] < 1e-14:
            return V[:, :j + 1], H[:j + 1, :j]
        V[:, j + 1] = w / H[j + 1, j]

    return V, H

大型問題一般只提供 $v\mapsto Av$ ，不形成完整矩陣。Restarted Arnoldi、shift-and-invert 與預條件可針對所需譜區域。

應用

結構工程

無阻尼自由振動

M\ddot x+Kx=0

導出廣義特徵值問題

Kx=\omega^2Mx.

$\omega$ 是自然頻率，特徵向量是模態形狀。相關方法亦用於屈曲臨界負載、氣動彈性 flutter 與地震反應。

量子化學

電子結構計算包含

H\psi=E\psi,

其中 $H$ 是 Hamiltonian， $E$ 是能階。以下只示意在指定基底建立矩陣：

def build_hamiltonian(basis_size, potential):
    H = np.zeros((basis_size, basis_size))
    for i in range(basis_size):
        H[i, i] = (i + 0.5) * np.pi**2 / 2
    for i in range(basis_size):
        for j in range(basis_size):
            H[i, j] += potential_matrix_element(i, j, potential)
    return H

它並非完整量子化學模型；基底、邊界條件、單位及 potential_matrix_element 必須由具體物理問題定義。

機器學習

PCA 對協方差矩陣作特徵分解，以主要特徵向量降維；譜聚類則利用圖 Laplacian 的低階特徵向量表示群落結構。高維情況通常使用 randomized SVD 或迭代方法，而非完整分解。

網絡分析

PageRank、eigenvector centrality、Katz centrality 與 hub-authority scores 都與譜結構有關。帶 teleportation 的 PageRank 可寫成線性系統

(I-\alpha P)^Tx=(1-\alpha)e,

但 $P$ 的行／列約定與 $e$ 的正規化必須保持一致。

大型問題的計算挑戰

大型譜問題需要：

記憶體管理
- matrix-free 乘法
- 分散式計算
- 階層式或稀疏儲存
效能
- 平行與 GPU 加速
- 預條件
- block Krylov 方法
- 通訊避免演算法

原文的 block_arnoldi 程式把 QR 的正交基錯誤存入 Hessenberg 陣列，維度亦不一致，因此不應當作可運行實作。正確 block Arnoldi 在每輪應先把 $W$ 對現有基底區塊正交化，再以 QR 得到新基底區塊 $V_{j+1}$ 與小矩陣區塊 $H_{j+1,j}$ ；實務上應使用 SciPy、SLEPc、ARPACK 或經驗證的專門函式庫。

未來方向

量子計算
- 量子相位估計及特徵值演算法
- 混合經典—量子方法
- 誤差緩解
機器學習整合
- 學習式預條件與演算法選擇
- 自適應精度
- 以代理模型降低昂貴矩陣運算
高效能計算
- GPU 與分散式方法
- fault-tolerant 實作
- 減少資料搬移與同步

結語

特徵值問題是數值分析的基石。由冪法走到帶位移 QR，再到 Lanczos 與 Arnoldi，方法演進反映同一原則：利用目標譜區域與矩陣結構，避免支付完整問題的成本。

理論界線、殘差、條件數及後向誤差與速度同樣重要。計算出很多小數位，不代表非正常或病態問題的特徵值已有同等精度。

參考文獻

Kincaid, D. & Cheney, W. (2009). Numerical Analysis: Mathematics of Scientific Computing, 3rd ed. AMS.
Golub, G. H. & Van Loan, C. F. (2013). Matrix Computations, 4th ed. Johns Hopkins University Press.
Trefethen, L. N. & Bau III, D. (1997). Numerical Linear Algebra. SIAM.
Parlett, B. N. (1998). The Symmetric Eigenvalue Problem. SIAM.
Stewart, G. W. (2001). Matrix Algorithms, Volume II: Eigensystems. SIAM.
Saad, Y. (2011). Numerical Methods for Large Eigenvalue Problems, 2nd ed. SIAM.