利用偽逆矩陣解線性方程

本文的閱讀等級:高級

在“線性代數的原罪?”一文裡,我曾說:

長久以來我們的數學教育方式 (其他課目也差不了多少) 是先告訴學生有關「數學」的事,展示給他們看這個「數學」是如何運作的,下課前塞給學生一些問題,回家自己去練習。

這種學習方式是學校教育講求成本效率下的產物,誰也說不準到底有多少長遠效果。我認為的理想學習方式應該帶有濃厚的實驗與實踐色彩,也就是經由發現問題,嘗試提出解答,從而建立理論。

 
舉例來說,有個學生正在進行一項研究工作,途中她面臨了線性聯立方程式的求解問題,A\mathbf{x}=\mathbf{b},如下:

\left[\!\!\begin{array}{rrr}    93.199&10.262&-29.334\\    1.560&32.957&61.377\\    27.152&37.982&57.482\\    23.434&-87.012&44.767    \end{array}\!\!\right]\begin{bmatrix}    x_1\\    x_2\\    x_3    \end{bmatrix}=\left[\!\!\begin{array}{r}    35.3000\\    73.4022\\    87.9935\\    -29.7854    \end{array}\!\!\right]

她察覺此方程式組總共有 4 條方程式,但只含 3 個變數,心生一計:「何不移除第四條方程式?這麼一來便可直接計算逆矩陣!」她決定將原方程式改為 \hat{A}\mathbf{x}=\hat{\mathbf{b}},如下:

\left[\!\!\begin{array}{rrr}    93.199&10.262&-29.334\\    1.560&32.957&61.377\\    27.152&37.982&57.482    \end{array}\!\!\right]\begin{bmatrix}    x_1\\    x_2\\    x_3    \end{bmatrix}=\begin{bmatrix}    35.3000\\    73.4022\\    87.9935    \end{bmatrix}

她又想到:「且慢,萬一這個 3\times 3 階係數矩陣 \hat{A} 是不可逆的,那不就糟了!」於是她計算了係數矩陣的行列式,得到 \mathrm{det}\hat{A}=-15.222\neq 0。確定了 \hat{A} 是可逆矩陣後,隨即利用本站提供的友好連結 Online Matrix Calculator 得到逆矩陣:

\hat{A}^{-1}=\left[\!\!\begin{array}{rrr}    28.694&111.946&-104.888\\    -103.589&-404.226&378.796\\    54.894&214.245&-200.732    \end{array}\!\!\right]

緊接著執行矩陣乘法運算得出解:

\mathbf{x}=\hat{A}^{-1}\hat{\mathbf{b}}=\hat{A}^{-1}\begin{bmatrix}    35.3000\\    73.4022\\    87.9935    \end{bmatrix}=\begin{bmatrix}    0.5129\\    0.8533\\    0.7247    \end{bmatrix}

為確定無誤,她又重新計算一次,這回貪圖方便僅輸入常數向量 \hat{\mathbf{b}} 至小數點後第二位,乘開發現

\mathbf{x}=\hat{A}^{-1}\hat{\mathbf{b}}^{\prime}=\hat{A}^{-1}\begin{bmatrix}    35.30\\    73.40\\    87.99    \end{bmatrix}=\begin{bmatrix}    0.6300\\    0.4302\\    0.9489    \end{bmatrix}

比較兩次運算結果,她嚇了一跳!為何兩次計算結果出現如此巨大的差異?她懷疑這不是軟體程式的運算誤差,而是常數向量捨去微小量所產生的錯誤。顯然錯誤又被係數矩陣放大,究竟 \hat{A} 有何神秘之處?

 
要回答這個問題,必須深入矩陣內部探究它的真實結構,我們採用的方法稱為“奇異值分解 (SVD)”。對 \hat{A} 執行 SVD,\hat{A}=\hat{U}\hat{\Sigma}\hat{V}^T,其中

\hat{\Sigma}=\begin{bmatrix}    100.012&0&0\\    0&99.947&0\\    0&0&0.002    \end{bmatrix}

由 SVD 可導出 \hat{A} 的逆矩陣

\hat{A}^{-1}=(\hat{U}\hat{\Sigma} \hat{V}^T)^{-1}=\hat{V}\hat{\Sigma}^{-1}\hat{U}^T

\hat{\Sigma}^{-1}=\begin{bmatrix}    0.01&0&0\\    0&0.01&0\\    0&0&500    \end{bmatrix}

問題就出在這裡,\hat{\Sigma}^{-1} 解釋了為何兩組解的差距 \Vert \hat{A}^{-1}\hat{\mathbf{b}}-\hat{A}^{-1}\hat{\mathbf{b}}^{\prime}\Vert 會這麼大。利用正交矩陣變換不改變向量長度的性質 (參見“特殊矩陣(3):么正矩陣(酉矩陣)”),就有

\Vert \hat{A}^{-1}(\hat{\mathbf{b}}-\hat{\mathbf{b}}^{\prime})\Vert=\Vert\hat{V}\hat{\Sigma}^{-1}\hat{U}^T(\hat{\mathbf{b}}-\hat{\mathbf{b}}^{\prime})\Vert=\Vert\hat{\Sigma}^{-1}\hat{U}^T(\hat{\mathbf{b}}-\hat{\mathbf{b}}^{\prime})\Vert

由於 \hat{A} 有一個數值很小的奇異值0.002\hat{U}^T(\hat{\mathbf{b}}-\hat{\mathbf{b}}^{\prime}) 的第三個元足足被放大 500 倍,這就是問題的根源!

 
要解決數值計算引發的不穩定問題,我們重新考慮原方程式組,計算原係數矩陣 A 的 SVD 分解 A=U\Sigma V^T,極小奇異值的情況不再發生:

\Sigma=\begin{bmatrix}    101.023&0&0\\    0&99.948&0\\    0&0&99.603\\    0&0&0    \end{bmatrix}

參考“通過推導偽逆矩陣認識線性代數的深層結構”,從 SVD 可求出 A 的偽逆矩陣:

A^{+}=V\Sigma^{+}U^{T}

其中 3\times 4\Sigma^{+} 的主對角元即為奇異值的倒數 1/\sigma_i,乘開後得到

A^{+}=\left[\!\!\begin{array}{rrrr}    0.009&0.000&0.003&0.002\\    0.001&0.003&0.004&-0.009\\    -0.003&0.006&0.006&0.005    \end{array}\!\!\right]

利用偽逆矩陣可算出方程式的解,如下:

\mathbf{x}=A^{+}\mathbf{b}= A^{+}\left[\!\!\begin{array}{r}    35.3000\\    73.4022\\    87.9935\\    -29.7854    \end{array}\!\!\right]=\begin{bmatrix}    0.5129\\    0.8533\\    0.7247    \end{bmatrix}

此解與先前 \hat{A}\mathbf{x}=\hat{\mathbf{b}} 的解相同,如果將常數向量 \mathbf{b} 的小數點後第三位元以下都捨去,結果與正確解相近:

\mathbf{x}= A^{+}\left[\!\!\begin{array}{r}    35.30\\    73.40\\    87.99\\    -29.79    \end{array}\!\!\right]=\begin{bmatrix}    0.5129\\    0.8533\\    0.7246    \end{bmatrix}

 
最後我們要問,給定一方陣 A 應該如何度量 \Vert A^{-1}\delta\mathbf{b}\Vert?其中 \delta\mathbf{b} 表示向量 \mathbf{b} 的微小變動量。這個問題的答案稱為條件數 (condition number),是數值線性代數裡的一個重要概念,將於另文詳細討論。

繼續閱讀:
Advertisements
This entry was posted in 線性代數專欄, 數值線性代數 and tagged , , , , . Bookmark the permalink.

4 則回應給 利用偽逆矩陣解線性方程

  1. levinc 說:

    這個例子再看一次還是覺得神奇~
    但有個地方用直觀理解卻怎也想不透,為什麼
    任意拿掉一條方程式計算的解,跟算完整的AX=b的解相同?

  2. ccjou 說:

    這是刻意挑選過的例子,四條方程式確實存在一解,也就是說,一條是多餘的。此例要說明的是 A^{+} 的數值穩定性比 A^{-1} 來的好,這反應在奇異值的大小。

  3. levinc 說:

    對,但那會是因為那最後一條方程式裡產生的問題嗎?能夠直觀的從原始資料矩陣結構觀察出嗎? 何不刪別條?

  4. ccjou 說:

    刪別條方程式可能不會發生病態系統
    https://ccjou.wordpress.com/2010/06/22/%E7%97%85%E6%85%8B%E7%B3%BB%E7%B5%B1/
    這由係數矩陣的奇異值決定,常用的檢查方法稱為條件數,從條件數大小即可判斷。
    https://ccjou.wordpress.com/2010/07/01/%E6%A2%9D%E4%BB%B6%E6%95%B8/

發表迴響

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

WordPress.com Logo

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

Twitter picture

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

Facebook照片

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

Google+ photo

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

連結到 %s