利用 Cayley-Hamilton 定理計算矩陣函數

本文的閱讀等級:中級

達文西說:「簡單是最終極的細緻。」(Simplicity is the ultimate sophistication.) 就數學而言,簡單不意味平庸,反而是優雅的體現。本文從一個簡單的問題開始,並致力於發展簡單的解法。令 A=\begin{bmatrix}  1&1\\  0&2  \end{bmatrix},求 A^{100}。典型的問題通常有典型的解法,對角化 (diagonalization) 是目前最常用的冪矩陣算法。矩陣 A 有相異特徵值 12,對應特徵向量 \begin{bmatrix}  1\\  0  \end{bmatrix}\begin{bmatrix}  1\\  1  \end{bmatrix},故可對角化為

A=\begin{bmatrix}  1&1\\  0&1  \end{bmatrix}\begin{bmatrix}  1&0\\  0&2  \end{bmatrix}\begin{bmatrix}  1&1\\  0&1  \end{bmatrix}^{-1}=S\Lambda S^{-1}

利用上式立得

A^{100}=(S\Lambda S^{-1})^{100}=S\Lambda^{100}S^{-1}=\begin{bmatrix}  1&1\\  0&1  \end{bmatrix}\begin{bmatrix}  1&0\\  0&2^{100}  \end{bmatrix}\begin{bmatrix}  1&1\\  0&1  \end{bmatrix}^{-1}=\begin{bmatrix}  1&2^{100}-1\\  0&2^{100}  \end{bmatrix}

不過,遺憾的是對角化並不適用於所有的矩陣。若 A=\begin{bmatrix}  1&2\\  0&1  \end{bmatrix},則 A 有兩個特徵值 1,但其特徵空間僅含一線性獨立向量 \begin{bmatrix}  1\\  0  \end{bmatrix},即知 A 不可對角化 (見“可對角化矩陣與缺陷矩陣的判定”)。針對不可對角化矩陣,典型的方法是分解出 A 的 Jordan 形式 A=MJM^{-1},接著計算 A^{100}=MJ^{100}M^{-1} (見“利用 Jordan form 解差分方程與微分方程”)。表面上,計算冪矩陣並不很困難,使用 Jordan 典型形式似乎有些「小題大作」。倘若不採用 Jordan 分解,那麼還有其他方法嗎?繼續閱讀前,建議讀者先花個幾分鐘想一想。

 
最簡單的辦法是修改 A 矩陣使之可對角化。我們可以在 A 的主對角元加入擾動量,譬如 A_{\epsilon}=\begin{bmatrix}  1&2\\  0&1+\epsilon  \end{bmatrix},其中 \epsilon 是一極小的數。如此一來,A_{\epsilon} 有相異特徵值 11+\epsilon,因此可對角化為

\displaystyle  A_{\epsilon}=\begin{bmatrix}  1&2\\  0&1+\epsilon  \end{bmatrix}=\begin{bmatrix}  1&2\\  0&\epsilon  \end{bmatrix}\begin{bmatrix}  1&0\\  0&1+\epsilon  \end{bmatrix}\frac{1}{\epsilon}\begin{bmatrix}  \epsilon &-2\\  0&1  \end{bmatrix}

按照先前方式計算冪矩陣,如下:

\displaystyle\begin{aligned}  A_{\epsilon}^{100}&=\begin{bmatrix}  1&2\\  0&\epsilon  \end{bmatrix}\begin{bmatrix}  1&0\\  0&(1+\epsilon)^{100}  \end{bmatrix}\frac{1}{\epsilon}\begin{bmatrix}  \epsilon &-2\\  0&1  \end{bmatrix}\\  &=\frac{1}{\epsilon}\begin{bmatrix}  1&2\\  0&\epsilon  \end{bmatrix}\begin{bmatrix}  1&0\\  0&1+100\epsilon+o(\epsilon^2)  \end{bmatrix}\begin{bmatrix}  \epsilon &-2\\  0&1  \end{bmatrix}\\  &=\frac{1}{\epsilon}\begin{bmatrix}  1&2+200\epsilon+o(\epsilon^2)\\  0&\epsilon+100\epsilon^2+o(\epsilon^3)  \end{bmatrix}\begin{bmatrix}  \epsilon &-2\\  0&1  \end{bmatrix}\\  &=\frac{1}{\epsilon}\begin{bmatrix}  \epsilon&200\epsilon+o(\epsilon^2)\\  0&\epsilon+100\epsilon^2+o(\epsilon^3)  \end{bmatrix}\\  &=\begin{bmatrix}  1&200+o(\epsilon)\\  0&1+100\epsilon+o(\epsilon^2)  \end{bmatrix},\end{aligned}

上式中 o(\epsilon^k) 代表最小次為 k 的多項式。令 \epsilon\to 0,則 A_{\epsilon}^{100}\to A^{100}=\begin{bmatrix}  1&200\\  0&1  \end{bmatrix}。這個方法的主要缺點在於引入代數符號運算。當 A 的尺寸增大時,計算過程變得尤其複雜。例如,A=\begin{bmatrix}  1&2&0\\  0&1&2\\  0&0&1  \end{bmatrix},令 \epsilon_1\epsilon_2 代表兩個相異極小數,不難想像 A_{\epsilon_1,\epsilon_2}=\begin{bmatrix}  1&2&0\\  0&1+\epsilon_1&2\\  0&0&1+\epsilon_2  \end{bmatrix} 有相當繁瑣的對角化展開式。

 
除了對角化,還有別的方法嗎?欲跳脫思路的陷阱,我們應暫時遺忘可對角化與不可對角化之間的對立,並將注意力集中在問題本身──冪矩陣。仔細觀察後發現 A=\begin{bmatrix}  1&2\\  0&1  \end{bmatrix} 滿足 (A-I)^2=0。這個事實稱為 Cayley-Hamilton 定理,它指出:對於任一 n\times n 階矩陣 A,令 p(t)=\det(A-tI)A 的特徵多項式 (或 p(t)=\det(tI-A) 亦可),則 p(A)=0,也就是說,矩陣 A 被其特徵多項式消滅 (見“Cayley-Hamilton 定理”)。此例 A 的特徵多項式是

p(t)=\begin{vmatrix}  1-t&2\\  0&1-t  \end{vmatrix}=(1-t)^2=t^2-2t+1

就有 p(A)=(A-I)^2=A^2-2A+I=0。利用恆等式 A^2=2A-I,我們可以計算更高次的冪矩陣:

\begin{aligned}  A^3&=A(A^2)=A(2A-I)=2A^2-A=2(2A-I)-A=3A-2I,\\  A^4&=A(A^3)=A(3A-2I)=3A^2-2A=3(2A-I)-2A=4A-3I,\\  \vdots  \end{aligned}

重複此過程,任一冪矩陣 A^k 皆可以表示成 AI 的線性組合。從另一個角度來看,設想 A^{100} 除以 (A-I)^2 的表達式:

A^{100}=q(A)(A-I)^2+c_1A+c_0I

其中 q(A) 代表商,c_1c_0 是未定係數。因為 (A-I)^2=0,可知 A^{100}=c_1A+c_0I。接下來要解出 c_1c_0。考慮多項式除法:

t^{100}=q(t)(t-1)^2+c_1t+c_0

為了從 t^{100} 的表達式得到兩個條件式 (因為有兩個未知係數),求上式等號兩邊的導數,即有

100t^{99}=q^{\prime}(t)(t-1)^2+2q(t)(t-1)+c_1=\left(q^{\prime}(t)(t-1)+2q(t)\right)(t-1)+c_1

t=1 代入上面兩式可得

c_1+c_0=1,~~c_1=100

解出 c_1=100c_0=-99,所以 A^{100}=100A-99I=\begin{bmatrix}  1&200\\  0&1  \end{bmatrix}

 
除了計算冪矩陣,Cayley-Hamilton 定理是否也可以用來計算一般的矩陣函數 f(A)?(關於矩陣函數的定義,請見“矩陣函數 (上)”和“矩陣函數 (下)”。) 答案是肯定的。利用泰勒展開式和 Cayley-Hamilton 定理 p(A)=0,任何 n\times n 階矩陣 (不論是否可對角化) A 的良好定義函數 f(A) 都可表示成矩陣多項式:

f(A)=q(A)p(A)+c_{n-1}A^{n-1}+\cdots+c_1A+c_0I=c_{n-1}A^{n-1}+\cdots+c_1A+c_0I

A 有相異特徵值 \lambda_1,\ldots,\lambda_m,且 \beta_1,\ldots,\beta_m 為對應的相重數 (亦稱代數重數),則 A 的特徵多項式為

p(t)=(t-\lambda_1)^{\beta_1}\cdots(t-\lambda_m)^{\beta_m}

考慮

\begin{aligned}  f(t)&=q(t)p(t)+c_{n-1}t^{n-1}+\cdots+c_1t+c_0\\  &=q(t)(t-\lambda_1)^{\beta_1}\cdots(t-\lambda_m)^{\beta_m}+c_{n-1}t^{n-1}+\cdots+c_1t+c_0,  \end{aligned}

對於特徵值 \lambda_i,我們可以使用 \beta_i 個限制條件:f(\lambda_i),~f^{\prime}(\lambda_i),\ldots,~f^{(\beta_i-1)}(\lambda_i)。因為 \sum_{i=1}^m\beta_i=n,所以共有 n 個條件式可以用來決定係數 c_0,c_1,\ldots,c_{n-1},以矩陣形式表示如下:

\begin{bmatrix}  1&\lambda_1&\lambda_1^2&\lambda_1^3&\cdots&\lambda_1^{n-1}\\  \vdots&\vdots&\vdots&\vdots&&\vdots\\  1&\lambda_m&\lambda_m^2&\lambda_m^3&\cdots&\lambda_m^{n-1}\\  -&-&-&-&-&-\\  \vdots&\vdots&\vdots&\vdots&&\vdots\\  0&1&2\lambda_i&3\lambda_i^2&\cdots&(n-1)\lambda_i^{n-2}\\  \vdots&\vdots&\vdots&\vdots&&\vdots\\  -&-&-&-&-&-\\  \vdots&\vdots&\vdots&\vdots&&\vdots\\  0&0&2&6\lambda_i&\cdots&(n-1)(n-2)\lambda_i^{n-3}\\  \vdots&\vdots&\vdots&\vdots&&\vdots\\  -&-&-&-&-&-\\  \vdots&\vdots&\vdots&\vdots&&\vdots  \end{bmatrix}\begin{bmatrix}  c_0\\  c_1\\  c_2\\  c_3\\  \vdots\\  \vdots\\  \vdots\\  c_{n-1}  \end{bmatrix}=\begin{bmatrix}  f(\lambda_1)\\  \vdots\\  f(\lambda_m)\\  -\\  \vdots\\  f^{\prime}(\lambda_i)\\  \vdots\\  -\\  \vdots\\  f^{\prime\prime}(\lambda_i)\\  \vdots\\  -\\  \vdots  \end{bmatrix}

其中由特徵值構成的 n\times n 階係數矩陣是可逆的,因此保證 c_0,c_1,\ldots,c_{n-1} 唯一存在,也就有唯一的 f(A)。係數矩陣為可逆矩陣的理由如下:(一) 最上面分塊包含相異特徵值的 m 個列屬於 Vandermonde 矩陣 V 的一部分,因此是線性獨立的 (見“特殊矩陣 (8):Vandermonde 矩陣”);(二) 隨後的各分塊皆可表示為 VD,其中 D 是可逆對角矩陣,因此也是線性獨立的;(三) 各分塊的領先行位置相異,如第一分塊起始於第一行,第二分塊起始於第二行,故各分塊彼此獨立。

 
下面以一個例子展示計算步驟。考慮 A=\begin{bmatrix}  1&2&0\\  0&1&0\\  0&0&3  \end{bmatrix},求 e^{A}

(1) 求出特徵多項式 p(t)=(t-1)^2(t-3),可知 A 有特徵值 1, 1, 3

(2) 令 e^A=c_2A^2+c_1A+c_0I。寫出函數 f(t)=e^t 的多項式除法表達式:

f(t)=e^{t}=q(t)(t-1)^2(t-3)+c_2t^2+c_1t+c_0

f(t) 的導數,

f^{\prime}(t)=e^t=g(t)(t-1)+2c_2t+c_1

其中 g(t) 是整理合併後的多項式。

(3) 將 t=1t=3 代入 f(t),並將 t=1 代入 f^{\prime}(t),可得限制條件式:

\begin{aligned}  e&=c_2+c_1+c_0\\  e^3&=9c_2+3c_1+c_0\\  e&=2c_2+c_1,  \end{aligned}

或表示為矩陣形式:

\begin{bmatrix}  1&1&1\\  1&3&9\\  0&1&2  \end{bmatrix}\begin{bmatrix}  c_0\\  c_1\\  c_2  \end{bmatrix}=\begin{bmatrix}  e\\  e^3\\  e  \end{bmatrix}

由上式解出

\displaystyle  c_0=\frac{1}{4}e^3-\frac{3}{4}e,~~c_1=-\frac{1}{2}e^3+\frac{5}{2}e,~~c_2=\frac{1}{4}e^3-\frac{3}{4}e

也就得到

\displaystyle\begin{aligned}  e^A&=\left(\frac{1}{4}e^3-\frac{3}{4}e\right)A^2+\left(-\frac{1}{2}e^3+\frac{5}{2}e\right)A+\left(\frac{1}{4}e^3-\frac{3}{4}e\right)I\\  &=\left(\frac{1}{4}e^3-\frac{3}{4}e\right)\begin{bmatrix}  1&4&0\\  0&1&0\\  0&0&9  \end{bmatrix}+\left(-\frac{1}{2}e^3+\frac{5}{2}e\right)\begin{bmatrix}  1&2&0\\  0&1&0\\  0&0&3  \end{bmatrix}+\left(\frac{1}{4}e^3-\frac{3}{4}e\right)\begin{bmatrix}  1&0&0\\  0&1&0\\  0&0&1  \end{bmatrix}\\  &=\begin{bmatrix}  e&2e&0\\  0&e&0\\  0&0&e^3  \end{bmatrix}.\end{aligned}

 
相較於對角化 f(A)=Sf(\Lambda)S^{-1} 或 Jordan 典型形式 f(A)=Mf(J)M^{-1},利用 Cayley-Hamilton 定理計算矩陣函數 f(A) 有三個明顯的優點:(一) 不須計算特徵向量或廣義特徵向量;(二) 不須計算特徵向量矩陣 SM 的逆矩陣;(三) 不須計算 Jordan 矩陣 J,自然也就免除 f(J) 的繁瑣算式 (見“矩陣函數 (下)”)。然而,此法必須使用一般矩陣乘法算出所有的冪矩陣 A^kk<n。提醒讀者,二個 n\times n 階矩陣相乘需要 O(n^3) 個乘法運算,因此乘法運算總量是 O(n^4)。換句話說,本文介紹的這個方法僅適用於小尺寸矩陣。

Advertisements
本篇發表於 特徵分析, 線性代數專欄 並標籤為 , , , 。將永久鏈結加入書籤。

發表迴響

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

WordPress.com Logo

您的留言將使用 WordPress.com 帳號。 登出 / 變更 )

Twitter picture

您的留言將使用 Twitter 帳號。 登出 / 變更 )

Facebook照片

您的留言將使用 Facebook 帳號。 登出 / 變更 )

Google+ photo

您的留言將使用 Google+ 帳號。 登出 / 變更 )

連結到 %s