漫步線性代數七——特殊矩陣和應用

本篇文章有兩個目标。第一是解釋實際問題中大型線性方程組 Ax=b 的一種解法，事實是，工程或經濟學中大型和現實的問題能夠引導我們更深入了解這些知識。但是有一個很重要應用卻不需要大量的準備工作。

另一個目标是說明系數矩陣具有的一些特殊性質，為了友善我們用同一個應用進行講解。大型矩陣幾乎總是有一個清晰的模式-對稱和很多零元素。因為一個稀疏矩陣包含的資訊遠小于 n2 個，是以計算應該更快。我們将觀察帶狀矩陣，看看集中在對角線附近是如何加快消元的，為此我們将看到一個特殊的三對角矩陣。

看方程(6)中的矩陣，它是通過将微分方程變化為矩陣方程得到的。這是對每個 x 求u(x)的連續問題，很明顯計算機不能解決它，是以它必須近似為一個離散的問題-我們保留更多的未知變量，結果的精度就越好，當然計算代價也就越高。作為一個簡單但仍然具有代表性的連續問題，我們選擇微分方程

−d2udx2=f(x)0≤x≤1(1)

這是關于位置函數 u(x) 的線性方程，可以在解中加上任何組合 C+Dx 依然滿足要求，因為 C+Dx 的二階導為零，不影響結果。對于兩個任意常數 C,D 的不确定性，通過在區間的兩端添加一個邊界條件就能夠移除：

u(0)=0,u(1)=0(2)

這個結果是一個兩點邊值問題，描述的不是瞬變而是穩态現象-例如一根棒的溫度分布，它的一端固定為 00C 并且熱源為 f(x) 。

記住，我們的目标是産生一個離散的問題-換句話說，一個線性代數中的問題。為此我們隻可以接受 f(x) 有限的資訊，描述在 n 個相等的區間點x=h,x−2h,…,x=nh上的值，對于同樣位置處的真是解 u 我們計算近似解u1,…,un，在端點處 x=0,x=1=(n+1)h 處，邊界值是 u0=0,un+1=0 。

第一個問題是：我們如何替換導數 d2u/dx2 ？一階導數可以近似表示為有限步長内停止的 Δu/Δx 并且不允許 h(orΔx) 趨近于零， Δu 可以是前面的，後面的或中間的：

ΔuΔx=u(x+h)−u(x)h or u(x)−u(x−h)h or u(x+h)−u(x−h)2h(3)

最後一個關于 x 對稱，它是最精确的。對于二階導數，隻是利用x,x±h處數值的一個組合：

d2udx2≈Δ2uΔx2=u(x+h)−2u(x)+u(x−h)h2(4)

它也有關于 x 對稱的優點。重複一遍：當h→∞時右邊接近 du/dx2 的真實值，但是我們必須讓 h 停在某個正數上。

對每個點x=jh，方程 −d2u/dx2=f(x) 可以用它的離散模拟(5)代替，我們通過乘以 h2 來得出 n 個方程Au=b：

−uj+1+2uj−uj−1=h2f(jh)for j=1,…,n(5)

第一個和最後一個( j=1,j=n )包含 u0=0,un+1=0 ，他們是已知的邊界條件，如果這些值非零的話，他們就轉化成右邊的值。這 n 個方程(5)的結構可以更矩陣形式來更好的可視化，我們選擇h=16，進而得到 5×5 的矩陣：

⎡⎣⎢⎢⎢⎢⎢⎢2−1−12−1−12−1−12−1−12⎤⎦⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢u1u2u3u4u5⎤⎦⎥⎥⎥⎥⎥⎥=h2⎡⎣⎢⎢⎢⎢⎢⎢f(h)f(2h)f(3h)f(4h)f(5h)⎤⎦⎥⎥⎥⎥⎥⎥(6)

現在，我們将求解方程(6)，它的系數矩陣非常有規律，有許多特殊的性質，其中有三個是非常基本的：

1. 矩陣 A 是三對角的。所有非零元素位于主對角線以及附近的兩條對角線上，這條窄帶以外的aij=0，這些零大大簡化了高斯消元過程。
2. 矩陣是對稱的。每個元素 aij 等于它的鏡像 aji ，使得 AT=A 。上三角矩陣 U 将是下三角矩陣L的轉置， A=LDLT 。 A 的對稱性反映了d2u/dx2的對稱性，奇導數像 du/dx,d3u/dx3 将破壞對稱性。
3. 矩陣是正定的。這個額外的性質說明主元是正的，在理論和實踐中不需要行變換。這和本文末尾要将的矩陣 B (非正定的)正好相反，在沒有行變換的情況下，它将對舍入非常敏感。另外關于正定的概念我會在以後的文章中詳細介紹！

我們傳回到A是三對角矩陣這個事實，它對消元會有什麼影響呢？消元過程的第一步是在第一個主元下面産生零：

⎡⎣⎢⎢⎢⎢⎢⎢2−1−12−1−12−1−12−1−12⎤⎦⎥⎥⎥⎥⎥⎥→⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢20−132−1−12−1−12−1−12⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥

跟一般的 5×5 矩陣相比，這一步主要有兩個簡化：

1. 在主元下面隻有一個非零元素
2. 主元所在的行非常短

乘數因子是 ℓ=−12 ，新的主元是 32 。更進一步，三對角模式保留着：每步消元都允許這兩個簡化。

最終結果是 LDU=LDLT ，注意主元！

A=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢1−121−231−341−451⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢2132435465⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢1−121−231−341−451⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

三對角矩陣的 L,U 分解因子是二對角矩陣，三個因子和矩陣 A 一樣對角線有同樣的帶狀結構，注意L,U互相之間存在轉置關系，我們從對稱就預期到會如此，主元 2/1,3/2,4/3,5/4,6/5 都是正的，他們的乘積就是 A 的行列式detA=6。當 n 不斷變大時，這些主元明顯收斂到1，這樣的矩陣計算起來非常友善。

稀疏因子L,U完全改變了通常的操作次數，每列的消元隻需要兩步，對于 n 列來說，次數從n3/3降到了 2n ，三對角方程組幾乎能很快解決，求解它的代價和 n 成正比。

帶狀矩陣就是在|i−j|<w内 aij=0 (圖1)，對角矩陣的半個帶寬 w=1 ，三對角矩陣的為 w=2 ,，當 w=n 時就是就是一般的矩陣了。對每一列，消元法需要 w(w−1) 次操作：長為 w 的行對下面的w−1行進行操作，對于 n 列的帶狀矩陣大約需要w2n次操作。

當 w 趨近n時，矩陣變成一般矩陣，次數變成 n3 。産生 L,D,U 的除法和乘-減(不考慮 A 為對稱這個假設)精确次數是P=13w(w−1)(3n−2w+1)，對于一般矩陣即 w=n ， P=13n(n−1)(n+1) ，這是一個整數，因為 n−1,n,n+1 是連貫的，是以他們之中有一個能被3整除。

圖1：帶狀矩陣 A 和它的因子L,U

這是我們最終得到的操作次數，我們強調一點，像 A 那樣的有限差分矩陣有逆，在求解Ax=b時，知道 A−1 比知道 L,U 要糟糕，因為 A−1 乘 b 需要n2步，而前向消元和回代(産生 x=U−1c=U−1L−1b=A−1b ) 4n 步就足夠了。

我們希望本例增強讀者對消元的了解，它是實踐中遇到的一個大型線性方程執行個體，接下來對于 n 個未知量的m個方程，我們将轉向讨論 x 的存在性和唯一性。

舍入誤差

理論上非奇異情況有完整的主元(考慮行變換)，實踐中，需要更多的行交換否則計算的解可能變得毫無價值。接下來我們重點來講如何使消元更加穩定-為什麼需要它以及如何做。

對于中等大小的系統，比如說100×100，消元可能涉及一百萬的三分之一次( 13n3 )，對于每步操作，我們必須預期舍入誤差。通常情況下，我們固定有效位數，然後将兩個大小不同的數相加給出一個誤差：

.456+.00123→.457丢掉數字2和3

那麼所有這些誤差是如何影響 Ax=b 的最終誤差的呢？

這不是一個簡單的問題。而這個問題被約翰 ⋅ 馮 ⋅ 諾依曼碰到了，在計算機使得100萬步操作成為可能的時候，是他引領了數學。事實上高斯和馮 ⋅ 諾依曼的結合給出了簡單的消元算法，雖然馮諾依曼高估了最後的舍入誤差。威爾金森(Wilkinson)找到這個問題的正确方法，并且他的書到現在依然是經典。

我們将給出兩個簡單的例子來說明舍入誤差的重要性：

A=[1.1.1.1.0001]B=[.0011.1.1.]

A 幾乎是奇異的而B離奇異很遠，如果我們稍微将 A 中的元素改一下a22=1，它既是奇異的。考慮兩個非常接近的向量 b ：

ill−conditioned uu++v1.0001v==22andwell−conditioned uu++v1.0001v==22.0001

第一個解是 u=2,v=0 ，而第二個解是 u=v=1 。 b 中第五位數的改變放大到解中第一位數的改變，沒有任何數值方法可以避免它對小擾動這種敏感，病态條件能夠從一個地方轉到另一個地方，但是無法移除。真實解都如此敏感更别說計算解了。

第二點如下：

15、即便是一個well−conditioned矩陣 B 也可能别差的算法個毀掉

我們很遺憾的說：對于矩陣B，直接使用高斯消元就是一種差的算法。假設.0001作為第一個主元，那麼第二行需要減去第一行的10000倍，右下角就變成了-9999，但是第三位四舍五入将變為-10000，為1的元素将會消失：

.0001u+vu+v==22→.0001u+v−9999v==1−9998

四舍五入将得到 −10000v=−10000 或者 v=1 。首先我們用正确值 v=.999 (保留三位小數)回代将得到 u=1 ：

.0001u+.9999=1oru=1

但是如果接受 v=1 ，我們将得到 u=0 ：

.0001u+1=1oru=0

計算結果完全不對， B 是well-conditioned但是消元明顯不穩定，L,D,U也明顯是比 B 大許多：

B=[11000001][.000100−9999][10100001]

小主元.0001帶來了不穩定，補救方法很明顯-換行

16、一個小主元迫使消元中發生變化。通常我們比較同一列所有可能的主元，将最大主元換到目前的位置

對于 B 來講，.0001将和下面的1比較，立馬進行行交換，從矩陣的角度來講就是乘以一個置換矩陣P，新矩陣 C=PB 就有好的因子：

C=[1.000111]=[1.000101][100.9999][1011]P=[0110]

C 的主元是1和.9999，比B的.0001和-9999要好。

還有一種政策就是與其餘所有列中最大主元進行交換，這時候可能不僅是行，也會有列交換。(這時候置換矩陣乘在右邊)這麼做的代價太高，上面的方法其實已經夠了。

我們說完了數值線性代數的基本算法：帶有行變換的消元法。一些進一步的完善，比如看整行或列，都是有可能的，但本質上讀者現在知道了一台電腦如何解線性方程組，相比“理論”描述-找到 A−1 并相乘 A−1b -我們的描述已花費了大量時間和耐心，我希望能夠用更簡單的方法來解釋 x <script type="math/tex" id="MathJax-Element-7260">x</script>是如何發現的，但是我認為我還沒有找到。(部落客心聲：期待大家能有更好的方法提出來)