Power 迭代法

Posted on 11/02/2010 by ccjou

本文的閱讀等級：中級

在線性代數的發展過程中，設計效率高且穩定性佳的矩陣特徵值算法是一個非常重要的研究問題。給定一個 n 階方陣 A，特徵值 \lambda 和對應的特徵向量 \mathbf{x}，\mathbf{x}\neq\mathbf{0}，滿足 A\mathbf{x}=\lambda\mathbf{x}，或寫為 (A-\lambda I)\mathbf{x}=\mathbf{0}，因此推論 \mathrm{det}(A-\lambda I)=0。表面上，欲決定矩陣 A 的特徵值 \lambda，只需要解出特徵多項式 p_A(t)=\mathrm{det}(A-tI) 的根即可。不幸的是，這個方法僅適用於小型矩陣，因為一般 n>4 次方程式並不存在公式解。針對高次多項式，雖然確實有可用的迭代 (iteration) 求根算法，然而，求多項式根在本質上是一個病態問題，也就是說，根的位置很容易受到多項式係數的微小擾動而發生劇烈改變。Wilkinson 多項式說明了這種情況，考慮

\begin{aligned} w(x)&=(x-1)(x-2)\cdots(x-20)\\ &=x^{20}-210x^{19}+20,615x^{18}-1,256,850x^{17}+53,327,946x^{16}-\cdots\end{aligned}

當 x^{19} 的係數從 -210 減小為 -210.0000001192，原本一根 x=20 改變為 x\approx 20.8，而 x=18 和 x=19 兩根則撞擊在一起成為重根 x\approx 18.62[1]。更有甚者，縱使原本是良置的特徵值問題也可能產生病態多項式，因為這個緣故，我們勢必要尋求其他數值計算方法。本文介紹一個簡易的特徵值算法，稱為 Power 迭代法，以及此法的一個改良版本。

 
假如方陣 A 絕對值最大的特徵值唯一存在， \vert\lambda_1\vert>\vert\lambda_2\vert\ge\cdots\ge\vert\lambda_n\vert\ge 0，Power 迭代法的演算程序如下：

Power 迭代法

  1. 給定一初始向量 \mathbf{u}_0
  2. 對於 k=0,1,2,\ldots，直至 \lambda(k) 收斂止，計算
  3.     \mathbf{u}_{k+1}=\frac{A\mathbf{u}_k}{\Vert A\mathbf{u}_k\Vert}
  4.     \lambda(k+1)=\mathbf{u}^T_{k+1}A\mathbf{u}_{k+1}

 
Power 迭代法的計算原理建立在差分方程上。當 A 是可對角化時，存在一組完整的線性獨立特徵向量 \mathbf{x}_1,\ldots,\mathbf{x}_n 對應特徵值 \lambda_1,\ldots,\lambda_n。若 A 是不可對角化矩陣，只要在 A 的各個元滲入相當捨入 (roundoff) 誤差的極小隨機擾動量，便能產生微小的特徵值變化，相同特徵值發生的機率也就幾近於零。因此就數值計算而言，所有的矩陣都是可對角化的。既然如此，初始向量可以唯一表示為 \mathbf{u}_0=c_1\mathbf{x}_1+c_2\mathbf{x}_2+\cdots+c_n\mathbf{x}_n，差分方程 \mathbf{u}_k=A^k\mathbf{u}_0 的解則為

\mathbf{u}_k=c_1\lambda_1^k\mathbf{x}_1+ c_2\lambda_2^k\mathbf{x}_2+\cdots+c_n\lambda_n^k\mathbf{x}_n

若 c_1\neq 0，亦即初始向量 \mathbf{u}_0 包含某些成分的 \mathbf{x}_1 (隨機選擇的 \mathbf{u}_0 通常滿足上述條件)，就有

\displaystyle\mathbf{u}_k=\lambda_1^k\left[c_1\mathbf{x}_1+c_2\left(\frac{\lambda_2}{\lambda_1}\right)^k\mathbf{x}_2+\cdots+ c_n\left(\frac{\lambda_n}{\lambda_1}\right)^k \mathbf{x}_n\right]

對於 j=2,\ldots,n，隨著 k 增大，\left(\lambda_j/\lambda_1\right)^k 逐漸逼近零，得知 \mathbf{u}_k\approx\lambda_1^kc_1\mathbf{x}_1。在不失一般性的情況下，假設所有的 \mathbf{x}_i 皆為單位向量 (unit vector)，我們可以推斷以下結果：

\displaystyle\mathbf{u}_{k+1}=\frac{A\mathbf{u}_k}{\Vert A\mathbf{u}_k\Vert}\approx\frac{\lambda_1^kc_1\mathbf{x}_1}{\Vert\lambda_1^kc_1\mathbf{x}_1\Vert}=\pm\mathbf{x}_1

\lambda(k+1)=\mathbf{u}^T_{k+1}A\mathbf{u}_{k+1}\approx(\pm\mathbf{x}_1)^TA(\pm\mathbf{x}_1)=\lambda_1

數列 \lambda(k) 將會收斂至絕對值最大的特徵值 \lambda_1，而向量序列 \mathbf{u}_k 則收斂至對應的特徵向量 \pm\mathbf{x}_1。

 
Power 法的優點是每個計算步驟僅需使用一次矩陣—向量乘法 \mathbf{y}=A\mathbf{u}_k；若將正規化特徵向量 \mathbf{u}_{k+1}=\mathbf{y}/\Vert\mathbf{y}\Vert 的分母改為 \Vert\mathbf{y}\Vert_{\infty}=\mathrm{max}_i\vert y_i\vert，還可以進一步減少計算量。Power 法的收斂速度由比值 \vert\lambda_2/\lambda_1\vert 決定，如果 \lambda_2 十分接近 \lambda_1，這將使得收斂緩慢。此外，Power 法最主要限制在於僅能算出最大特徵值和對應的特徵向量，若要得到其他特徵值和特徵向量，此法還必須加以修改。

 
Power 迭代法的改良版稱為逆迭代法 (inverse iteration method)，或移動逆迭代法 (shift inverse method)，它與 Power 法的主要不同處在於以 (A-\alpha I)^{-1} 替換 A，詳細演算程序如下：

Power 移動逆迭代法

  1. 給定一初始向量 \mathbf{u}_0
  2. 對於 k=0,1,2,\ldots，直至 \lambda(k) 收斂止，計算
  3.     \mathbf{y}=(A-\alpha I)^{-1}\mathbf{u}_k
  4.     \mathbf{u}_{k+1}=\frac{\mathbf{y}}{\Vert\mathbf{y}\Vert}
  5.     \lambda(k+1)=\mathbf{u}^T_{k+1}A\mathbf{u}_{k+1}

 
純量 \alpha 稱為移動量 (shift)，此值決定上述程序最後收斂至那一個特徵值。假設 \alpha 和所有的特徵值都不相等，而且 \alpha 與特徵值 \lambda_j 的距離最近，可知 0<\vert\lambda_j-\alpha\vert<\vert\lambda_i-\alpha\vert，i\neq j。由 A\mathbf{x}_i=\lambda_i\mathbf{x}_i 可推得 (A-\alpha I)\mathbf{x}_i=(\lambda_i-\alpha)\mathbf{x}_i，\lambda_i-\alpha 為 A-\alpha I 的特徵值。因為 \alpha\neq\lambda_i，A-\alpha I 是可逆矩陣，故 (A-\alpha I)^{-1} 有特徵值 1/(\lambda_i-\alpha)，對應的特徵向量仍為 \mathbf{x}_i。根據 Power 法的結論，當 k 逐漸增大，\lambda(k) 將收斂至絕對值最大的特徵值，亦即 1/(\lambda_j-\alpha)，而 \mathbf{u}_k 則收斂至特徵向量 \mathbf{x}_j。

 
我們可以由 Gershgorin 圓估計出的特徵值範圍設定鄰近 \lambda_j 的一個固定 \alpha 值 (見“估計特徵值範圍的 Gershgorin 圓”)，另一個方法是令 \alpha 也隨迭代程序靠近 \lambda_j，如下：

\alpha=\lambda(k+1)= \mathbf{u}^T_{k+1}A\mathbf{u}_{k+1}

還有，計算 \mathbf{y}=(A-\alpha I)^{-1}\mathbf{u}_{k} 並不需要真的得到逆矩陣，解線性方程 (A-\alpha I)\mathbf{y}=\mathbf{u}_{k} 即可：計算 A-\alpha I 的 LU 分解得到 A-\alpha I=LU，以正向替代解 L\mathbf{z}=\mathbf{u}_k，之後再以反向替代解 U\mathbf{y}=\mathbf{z}。

 
Power 法和逆迭代法的主要應用在於計算大型稀疏矩陣 (包含許多零元) 的特徵值和特徵向量。如果我們想要得到完整的特徵值和特徵向量，或矩陣僅含少數的零元，這時便不宜採用 Power 迭代法。在一般情況下，QR 演算法是目前最常被採用的特徵值算法，這將留待日後再詳細介紹。

 
引用來源：
[1] http://en.wikipedia.org/wiki/Wilkinson%27s_polynomial

This entry was posted in 線性代數專欄, 數值線性代數 and tagged , , . Bookmark the permalink.

2 Responses to Power 迭代法

  1. Caren Guo says:
    11/03/2015 at 1:37 pm

    What if the matrix A has several same eigenvalues (for example, A has a Jordan block of (\lambda, 1), (0, \lambda))? What will the convergence rate be?

    Reply
    • ccjou says:
      11/03/2015 at 3:16 pm

      Let A=\begin{bmatrix} \lambda&1\\ 0&\lambda \end{bmatrix} and u_0=c_1e_1+c_2e_2, where e_1 and e_2 are standard unit vectors. If c_2=0, the process converges immediately. Consider the case of c_2\neq 0. Then,
      u_k=A^ku_0=\begin{bmatrix} \lambda^k&k\lambda^{k-1}\\ 0&\lambda^k \end{bmatrix}(c_1e_1+c_2e_2) =k\lambda^{k-1}\left((c_1\lambda/k+c_2)e_1+(c_2\lambda/k)e_2\right).
      As k\to\infty, u_{k+1}=\frac{Au_k}{\Vert Au_k\Vert}\approx e_1. The sequence converges sublinearly.

      Note that the power method fails to converge if there are two different dominant eigenvalues with equal absolute value, for example, B=\begin{bmatrix} \lambda&0\\ 0&-\lambda \end{bmatrix}.

      Reply

Leave a comment