特徵值的擾動分析

本文的閱讀等級:高級

若論方陣最精華的本質,誠然非特徵值莫屬。特徵值決定了矩陣所代表的線性變換的固有特性。特徵值的絕對值 (或稱向徑) 表示線性變換的伸縮大小,特徵值的幅角則表示旋轉強弱 (見“解讀複特徵值”)。若 A 的所有特徵值的絕對值小於 1,則冪矩陣 A^k 將隨 k 增大而收斂至零矩陣 (見“收斂矩陣”)。若 A 的所有特徵值的實部是負數,當 t\to\infty,矩陣指數 e^{At} 趨於零矩陣 (見“線性微分方程的穩定性”)。透過對矩陣特徵值的研究,不僅可以幫助我們了解矩陣的本質,還可以提供解讀複雜動態系統行為的線索。本文探討特徵值的擾動分析,也就是在引入擾動的情況下 (如數值計算的捨入誤差或來自線性系統外部的干擾),設法界定矩陣特徵值的變化範圍。我們將運用矩陣範數、對角化和三角化推導出兩個特徵值變異的上界[1]

 
Bauer-Fike 上界

A 是一 n\times n 階可對角化矩陣,A=S\Lambda S^{-1},其中 \Lambda=\hbox{diag}(\lambda_1,\ldots,\lambda_n)S 是特徵向量構成的可逆矩陣。為了方便,我們用 \sigma(A) 表示 A 的特徵值 \{\lambda_i\} 形成的集合,稱為矩陣譜 (spectrum)。令 B=A+E\mu\in\sigma(B),這裡 E 代表微小的擾動矩陣。下面結果稱為 Bauer-Fike 上界:

\displaystyle  \min_{\lambda_i\in\sigma(A)}\vert\mu-\lambda_i\vert\le\kappa_p(S)\Vert E\Vert_p

其中 \kappa_p(S)=\Vert S\Vert_p\cdot\Vert S^{-1}\Vert_p 稱為條件數 (condition number,見“條件數”)。上式中,\Vert\cdot\Vert_p 表示 p-範數,p 是 1, 2 或 \infty (見“矩陣範數”)。往後我們採用簡寫符號,除非特別聲明,\Vert\cdot\Vert 可為任一 p-範數。

 
如果 \mu\in\sigma(A),命題立刻成立。假設 \mu\notin\sigma(A),即知 (\mu I-A) 是可逆矩陣。考慮

(\mu I-A)^{-1}(\mu I-B)=(\mu I-A)^{-1}(\mu I-A-E)=I-(\mu I-A)^{-1}E

因為 \mu I-B 不可逆,故 I-(\mu I-A)^{-1}E 不可逆。利用 Neumann 無窮級數性質 (見“Neumann 無窮級數”):若 \Vert M\Vert<1,則 I-M 是可逆矩陣。我們斷定 \Vert(\mu I-A)^{-1}E\Vert\ge 1,否則 I-(\mu I-A)^{-1}E 的逆矩陣存在 。利用矩陣範數不等性質,\Vert XY\Vert\le\Vert X\Vert\cdot\Vert Y\Vert,可得

\displaystyle\begin{aligned}  1\le\Vert(\mu I-A)^{-1}E\Vert&=\Vert S(\mu I-\Lambda)^{-1}S^{-1}E\Vert\\  &\le\Vert S\Vert\cdot\Vert(\mu I-\Lambda)^{-1}\Vert\cdot\Vert S^{-1}\Vert\cdot\Vert E\Vert\\  &=\kappa(S)\Vert E\Vert\max_{\lambda_i\in\sigma(A)}\vert \mu-\lambda_i\vert^{-1}\\  &=\kappa(S)\Vert E\Vert\frac{1}{\min_{\lambda_i\in\sigma(A)}\vert\mu-\lambda_i\vert},\end{aligned}

此即所求。補充解釋 \Vert(\mu I-\Lambda)^{-1}\Vert=\max_{\lambda_i\in\sigma(A)}\vert\mu-\lambda_i\vert^{-1} 的由來:因為 (\mu I-\Lambda)^{-1} 是對角矩陣,不論何種 p-範數都等於最大的主對角元絕對值。

 
\kappa(S) 是一個相對小的數,則 A 的特徵值 \lambda_i 較不易受到微擾而變動,也就是相對不敏感。但如果 \kappa(S) 相對地大,我們就必須抱持謹慎的態度。若 A 是正規矩陣 (normal matrix,見“特殊矩陣 (2):正規矩陣”),則 S 是么正 (unitary) 矩陣 S^\ast=S^{-1},或正交 (orthogonal) 矩陣 S^T=S^{-1},就有 \Vert S\Vert=\Vert S^{-1}\Vert=1,故 \kappa(S)=1。在此情形下,特徵值變化的上界完全由微擾矩陣決定:\displaystyle  \min_{\lambda_i\in\sigma(A)}\vert\mu-\lambda_i\vert\le\Vert E\Vert,所以特徵值擾動分析通常僅針對非正規矩陣實施。

 
Schur 上界

使用 Schur 定理將矩陣 A 三角化為 A=U(\Lambda+N)U^\ast,其中 N=[n_{ij}] 是上三角矩陣且主對角元皆為零,U 是么正矩陣 (見“矩陣三角化的 Schur 定理”)。令 m 是最小正整數使得 \vert N\vert^m=0,這裡 \vert N\vert 表示絕對值矩陣,即 (\vert N\vert)_{ij}=\vert n_{ij}\vert。Schur 上界表達如下:

\displaystyle  \min_{\lambda_i\in\sigma(A)}\vert\mu-\lambda_i\vert\le\max\{\theta,\theta^{1/m}\}

其中 \theta=\Vert E\Vert_2\sum_{k=0}^{m-1}\Vert N\Vert_2^k

 
我們定義

\displaystyle  \delta=\min_{\lambda_i\in\sigma(A)}\vert\mu-\lambda_i\vert=\frac{1}{\Vert(\mu I-\Lambda)^{-1}\Vert_2}

同樣考慮 \mu\notin\sigma(A),則 \delta>0I-(\mu I-A)^{-1}E 不可逆。類似先前的推導方式,

\displaystyle\begin{aligned}  1\le\Vert(\mu I-A)^{-1}E\Vert_2&\le\Vert (\mu I-A)^{-1}\Vert_2\cdot\Vert E\Vert_2\\  &=\Vert ((\mu I-\Lambda)-N)^{-1}\Vert_2\cdot \Vert E\Vert_2.  \end{aligned}

這個結果將留待稍後使用。因為 (\mu I-\Lambda)^{-1} 是對角矩陣且 \vert N\vert^m=0,不難證明 ((\mu I-\Lambda)^{-1}N)^m=0。接下來我們運用一些矩陣代數技巧化簡。設 X 是可逆矩陣,考慮 X-Y=X(I-X^{-1}Y),即有 (X-Y)^{-1}=(I-X^{-1}Y)^{-1}X^{-1}。令 X=\mu I-\LambdaY=N。套用等式 1-a^m=(1-a)(1+a+\cdots+a^{m-1}),則有

\displaystyle  I=I-((\mu I-\Lambda)^{-1}N)^m=(I-(\mu I-\Lambda)^{-1}N)\sum_{k=0}^{m-1}((\mu I-\Lambda)^{-1}N)^k

故可得

\displaystyle\begin{aligned}  ((\mu I-\Lambda)-N)^{-1}&=(I-(\mu I-\Lambda)^{-1}N)^{-1}(\mu I-\Lambda)^{-1}\\  &=\sum_{k=0}^{m-1}((\mu I-\Lambda)^{-1}N)^k(\mu I-\Lambda)^{-1}.\end{aligned}

再用 \delta 的定義式和矩陣範數不等式,

\displaystyle  \Vert ((\mu I-\Lambda)-N)^{-1}\Vert_2\le\frac{1}{\delta}\sum_{k=0}^{m-1}\left(\frac{\Vert N\Vert_2}{\delta}\right)^k

分開兩種情況討論。若 \delta>1,則

\displaystyle  \Vert ((\mu I-\Lambda)-N)^{-1}\Vert_2\le\frac{1}{\delta}\sum_{k=0}^{m-1}\Vert N\Vert_2^k

由先前得到的不等式,可知

\displaystyle  \delta\le\Vert ((\mu I-\Lambda)-N)^{-1}\Vert_2^{-1}\sum_{k=0}^{m-1}\Vert N\Vert_2^k\le \Vert E\Vert_2\sum_{k=0}^{m-1}\Vert N\Vert_2^k=\theta

\delta\le 1,則

\displaystyle  \Vert ((\mu I-\Lambda)-N)^{-1}\Vert_2\le\frac{1}{\delta^m}\sum_{k=0}^{m-1}\Vert N\Vert_2^k

以相同方法可推得 \delta^m\le\theta。合併以上結果,\delta\le\max\{\theta,\theta^{1/m}\}

 
我們用一個例子比較兩種方法得到的上界 (取自[1]):

A=\begin{bmatrix}  1&2&3\\  0&4&5\\  0&0&4.001  \end{bmatrix},~~~E=\begin{bmatrix}  0&0&0\\  0&0&0\\  0.001&0&0  \end{bmatrix}

計算得到 A+E 的特徵值:\sigma(A+E)\approx\{1.0001, 4.0582, 3.9427\}。矩陣 A 的特徵向量矩陣 S 的條件數為 \kappa_2(S)\approx 10^7。Bauer-Fike 上界的數量級約為 10^4,但 Schur 上界的數量級約為 10^0。這兩個上界的差異在於 \kappa_2(S)\Vert N\Vert_2^{m-1},如果它們的數值偏大,縱使微小的擾動也可能引發特徵值的巨大改變。見下例:

A=\begin{bmatrix}  0&I_9\\  0&0  \end{bmatrix},~~~E=\begin{bmatrix}  0&0\\  10^{-10}&0  \end{bmatrix}

其中 E 的右上零分塊是 9\times 9 階。對於 \lambda\in\sigma(A)=\{0\} 以及所有 \mu\in\sigma(A+E)\vert \lambda-\mu\vert=\vert\mu\vert\approx 10^{-1}。這個例子顯示數量級為 10^{-10} 的一個微擾元就能造成數量級為 10^{-1} 的特徵值改變量。

 
參考來源:
[1] Gene H. Golub 和 Charles F. Van Loan 所著 Matrix Computations,第二版,1989年,頁341-345。

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

發表迴響

在下方填入你的資料或按右方圖示以社群網站登入:

WordPress.com Logo

你正使用 WordPress.com 帳號留言。 登出 / 變更 )

Twitter picture

你正使用 Twitter 帳號留言。 登出 / 變更 )

Facebook照片

你正使用 Facebook 帳號留言。 登出 / 變更 )

Google+ photo

你正使用 Google+ 帳號留言。 登出 / 變更 )

連結到 %s