1、第五章 矩陣代數(shù)數(shù)值計(jì)算,一、矩陣的基本運(yùn)算 二、矩陣的三角分解 三、矩陣的正交變換 四、矩陣的譜分解 五、IMSL中的線性系統(tǒng)、特征 值分析模塊,矩陣代數(shù)運(yùn)算是統(tǒng)計(jì)模型的基礎(chǔ)，統(tǒng)計(jì)模型的所有估計(jì)幾乎都是用矩陣代數(shù)運(yùn)算計(jì)算出結(jié)果。例如最小二乘估計(jì)、典型相關(guān)分析、因子分析以及各類回歸分析。從計(jì)算的角度來說為使計(jì)算結(jié)果可靠，我們總是先對矩陣進(jìn)行三角分解，然后進(jìn)行各種計(jì)算例如，矩陣的逆、求解線性方程組以及對矩陣進(jìn)行譜分解等。 本章首先介紹矩陣的三角分解，然后引導(dǎo)學(xué)習(xí)者使用IMSL和SASD中的豐富矩陣的算法，將它們拼接起來就可以解決各種矩陣的計(jì)算。,5.1 引言,矩陣代數(shù)運(yùn)算在數(shù)值計(jì)算中起著基礎(chǔ)性的

2、作用，無論我們建立了多么復(fù)雜的數(shù)學(xué)模型，最終我們總是要把它變?yōu)榫仃嚧鷶?shù)的形式。特別是統(tǒng)計(jì)模型，無論是多元線性回歸、廣義線性模型、多元統(tǒng)計(jì)分析無不與矩陣代數(shù)有著密切的聯(lián)系。我們所研究的對象，即樣本可以看成是一個(gè)矩陣，而統(tǒng)計(jì)上的協(xié)方差，實(shí)際上是該矩陣的轉(zhuǎn)置與該矩陣相乘形成的新矩陣的元素。而回歸的最小二乘估計(jì)的算法為：,（5.1.1),包含了矩陣的乘、矩陣的求逆及矩陣與向量的乘法等等。而特征值與特征向量在數(shù)理統(tǒng)計(jì)理都有明確的條件統(tǒng)計(jì)含義。因此我們將在這一章介紹矩陣運(yùn)算的基本數(shù)值方法，以及如何調(diào)用IMSL和SASD中豐富的算法模塊。,5.2 矩陣數(shù)值計(jì)算基礎(chǔ),對于一般的二維樣本，我們總可以寫成如下的矩

3、陣形式。,（5.2.1),從計(jì)算的角度處理矩陣問題的一個(gè)最有效的方法是，將一個(gè)一般的矩陣分解為幾個(gè)簡單矩陣的乘積形式。其中最便于計(jì)算的是三角矩陣，以上三角矩陣為例，由其性質(zhì)，兩個(gè)上三角陣的和、積仍為上三角陣，上三角陣的特征值就是其對角線元素。,系數(shù)矩陣為上三角陣的線性線性方程組是最容易求解的，上三角陣的逆陣仍然是上三角陣。因此處理矩陣計(jì)算問題的關(guān)鍵是將一般矩陣化為三角陣和對角陣的形式，然后進(jìn)行計(jì)算。 5.2.1 矩陣的三角 三角分解 （1）L*R分解,（5.2.2),其中 L 單位下三角陣（主對角線元素為1），R為上三角陣。 （2）LDR* 分解,(5.2.3),其中L為單位下三角陣，D為對角

4、陣，為單位上三角陣。,（3）Crout 分解,（5.2.4),其中為單位下三角陣，為單位上三角陣。 (4) Cholesky 分解,（5.2.5),其中A為正定對稱陣，T 為上三角陣。Cholesky 分解是統(tǒng)計(jì)計(jì)算中最常用的分解方法之一。因?yàn)槲覀兊膮f(xié)方差矩陣、相關(guān)矩陣都是使用這種分解方法。,5.2.2 矩陣的三角分解算法 以上四種分解是類似的，使用待定系數(shù)法。 (1) 以 LR* 分解為例，設(shè),其中A為正定陣，并記分解已為以下形式：,=,利用矩陣相對應(yīng)元素相等的事實(shí)，我們立即得到,現(xiàn)在我們可以計(jì)算矩陣 L 的第一列,由第二行，第二列相等,以及用前面的計(jì)算結(jié)果我們有：,矩陣 R* 的第二行,矩

5、陣L的第二列,即有以下公式：,從而我們可以推出一般的計(jì)算公式：,（2）Cholesky分解算法 同樣，利用待定系數(shù)法以及矩陣 A 的正定對稱性，我們有：,我們可以推導(dǎo)出 Cholesky 三角分解得算法 ：,為保證除法運(yùn)算時(shí) ，我們由以下定理 定理5.2.1 當(dāng) A 為對稱正定陣時(shí)，A 的Cholesky 分解必存在，并且當(dāng)限定 T 的對角元素為正時(shí)，其分解是唯一的。 有了矩陣三角三角分解后，各種矩陣的求解就十分方便了。例如：求解線性方程組,對 A 作 LR 分解，有 ，則解方程變換為解,5.2.3 矩陣的正交變換 我們從另一個(gè)角度來考慮LR分解，由前面的結(jié)論我們有,此表達(dá)式可以了解為對A線性

6、變換后變成了三角陣R，其中為變換陣。問題是我們能否用更為簡單的一系列變換將A變?yōu)樯先顷嚒?（1） 矩陣的正交三角分解 矩陣的正交分解可以寫成以下形式,其中 Q 是正交矩陣， 即 ，R是上三角矩陣，從而我們有,(5.2.6),這種變換在矩陣的運(yùn)算中是非常重要的。以下我們將 分解為一系列較為簡單的正交變換。 （2）Householder 變換 為產(chǎn)生盡量簡單的正交變換，我們考慮以下形式的正交方陣,(5.2.7),這里 In 是單位矩陣，u 為 n 維向量，為正實(shí)數(shù)。具有這種形式的正交變換稱為 H 型變換，我們可以通過以下步驟將矩陣 A變換為上三角陣 R，先用 H 型變換將A的第一列變量變?yōu)椋?再

7、用 H 型變換將A的第二列變?yōu)?：,第 i 步有,為實(shí)現(xiàn)這一過程，我們先考慮以下簡單問題。設(shè),我們要求一個(gè)H 型正交矩陣Hi ，使得后 n-I 個(gè)元素為 0，,其中 為常數(shù)， 為使后n-i個(gè)元素為0，可以取,這里,稱此 為由向量 定義的Householder變換，并有性質(zhì)。,1）,2） Hi 是正交的,3）,（3） Gives變換,Gives 變換具有以下性質(zhì)：,第 j 列,第 i 列,Gives 變換具有以下性質(zhì)： 1） 是正交矩陣,2）用 左 乘，結(jié)果只改變 A 的第 i 行第 j 行元素。,用 左 乘，結(jié)果只改變 A 的第 i 行第 j 行元素。,3）,對向量,這里,5.2.4 矩陣的譜

8、分解 前面的方法是用正交變換方法將矩陣A變?yōu)槿顷?，以下我們用同樣的方法將A變換為對角矩陣。 （1）對稱矩陣的譜分解,設(shè) 是 n 階方陣，以下分解式稱為 A 的譜分解式，或稱為 特征值分解式：,（5.2.8),其中 U 為 n 階正交方陣，D為對角陣， 稱 為矩陣 A 的譜或稱特征值，若記,記 ，則上式可以寫成,(5.2.9),如果 A 是實(shí)對稱矩陣，則 A 的譜分解一定存在。 （2）矩陣譜分解的計(jì)算方法,（5.2.9) 可以改寫如下：,（5.2.10),即A是經(jīng)過正交變換后化為對角陣的，我們可以利用Householder和Givens方法的思路來構(gòu)造這樣的正交變換，具體來講，我們可以將（5.

9、2.8）式中的U分解為一系列簡單的正交矩陣乘積的形式，具體算法為：,即,以下介紹如何適當(dāng)選取 ，使 在k充分大時(shí)接近于 一個(gè)對角陣。 Jacobi旋轉(zhuǎn)法,在 Gives 矩陣中取,其中 待定，對于 任意，可以驗(yàn)證 是一個(gè)正交變換，與解析幾何中的旋轉(zhuǎn)變換類似，在 n 維空間中，若對其中的二維(i, j)作旋轉(zhuǎn)變換，稱其為(i, j)平面上的旋轉(zhuǎn)矩陣。 可以證明 Jacobi 算法必有 為對角矩陣。 在（5.2.9）如果取 為 的 正交三角變換，則,為著名的求特征值的QR算法,3）QR算法 （a）,（b）對 做 的正交分解,取,A為任意方陣，可以證明 “基本收斂”于一個(gè)上三角矩陣， 而對角線元素為

10、其特征值。,5.2.5 矩陣的奇異值分解,設(shè) 是任意非零矩陣，則,為 A的奇異值分解，其中 U 和 V 分別為 n 和 m 階正交方陣， D 為 nm 階對角陣，其非對角元素均為 0，D的對角線元素 稱為矩陣 A 的奇異值。奇異值的分解與矩陣的譜分解方法類 似。,5.2.6 矩陣的廣義逆,矩陣的廣義逆在統(tǒng)計(jì)計(jì)算中具有重要的作用，它是由矩陣的逆 的概念進(jìn)一步一般化而來。 設(shè) 為 nm 階矩陣，G為 nm 矩陣，如果G滿足： （1）AGA=A （2）GAG=G （3）（AG）=AG（即AG為對稱矩陣） （4）（GA）=GA （即GA為對稱矩陣） 滿足這四個(gè)條件的某幾個(gè)或全部，則稱矩陣G為矩陣A的廣

11、義逆。,定義： 1）滿足上面第一條的矩陣 G 稱為 A的減號逆，記為 G=A- 2）滿足上面條件（1），（2）的矩陣 G 稱為 A的自反廣義 逆記為 3）滿足上面所有的四個(gè)條件稱G為A的加號逆極為,5.2.7 線性方程組的解 我們可以用矩陣的廣義逆這一有力的工具來解線性方程組,其中 A 為 nm 矩陣，b 為 n1的向量，x 為 m1的位知 向量，若（5.2.11）有解，則稱其為相容方程。,（5.2.11),如果存在 使得,則稱 為方程（5.2.11) 的最小二乘解，其中 表示向量 y 的模，其定義為。,以下不證明，給出相容方程的一般表達(dá)式,或,u任意,5.2.8 矩陣的范數(shù)（模）,矩陣的范數(shù)

12、與條件數(shù)是矩陣代數(shù)運(yùn)算的重要概念之一，范數(shù)為任意矩陣定義了一個(gè)函數(shù)，而條件數(shù)是將來進(jìn)行計(jì)算時(shí)對計(jì)算精度的一種衡量。請見范數(shù)的定義： 定義：設(shè) Rm 為 m 維實(shí)數(shù)空間， 是 Rm 到實(shí)數(shù)軸R1的一個(gè)映射，若 滿足 （1） 時(shí)成立 （2）對任意常數(shù) ， （3）,則稱 為向量 x 的范數(shù),常用的范數(shù)有： 1） 2） 3）,矩陣范數(shù)的定義： 定義：設(shè) 是nm維實(shí)數(shù)空間， 是 到實(shí)數(shù) 軸的一個(gè)映射，如果滿足：,（2）對任意常數(shù) ， 有,（3）,則稱 為矩陣A的范數(shù) ，特別，如果,（1）,（4）對,則稱 為矩陣A的相容范數(shù)。,矩陣A的常用范數(shù) 對于任意矩陣 定義：,容易證明 是矩陣A的相容范數(shù)。則常用的

13、矩陣范數(shù)：,（1）,（2）,（3）,這里我們主要討論 ，并記以,定理：設(shè) A 是 nm 矩陣，則有,這里 為矩陣 A 的絕對值最大的特征值（奇異值）， 從而我們給矩陣定義了一個(gè)范數(shù)。,5.2.9 矩陣的條件數(shù) 定義：稱 為 的條件數(shù)，記為：,條件數(shù)有以下性質(zhì)：,（1） 即矩陣的條件數(shù)為矩陣 的最大奇異值比最小奇異值的絕對值。特別當(dāng)A為對稱陣時(shí) 即為A的最大特征值和最小特征值之比的絕對值。,（2） ， 當(dāng) A 為正交陣時(shí)等式成立。,（3）對任意非零常數(shù)有，,（4）若 ，則,（5）,矩陣條件數(shù)的重要意義在于，可以判別一個(gè)求逆矩陣的病 態(tài)性。當(dāng)系數(shù)矩陣的條件數(shù)非常大時(shí)，即存在接近與0的特征 值，將導(dǎo)

14、致解的誤差急劇放大?；蛘哒f得出的解不可信。 對于最小二乘估計(jì)：,我們最關(guān)心的是對稱矩陣 的條件數(shù)是否正常，從而判 別最小二乘估計(jì)是否可信。,5.3 IMSL庫中的矩陣計(jì)算模塊,IMSL庫中的maths.lib庫中有非常豐富的矩陣計(jì)算模塊包括： 矩陣乘、矩陣求逆、廣義逆矩陣與向量乘等。 矩陣的譜分解、矩陣的奇異值分解。 矩陣、向量的模計(jì)算。 矩陣的條件數(shù)的計(jì)算。 求解線性方程組。,打開maths.lib我們可以看到,基本矩陣向量運(yùn)算模塊,求解線性方程組,矩陣的特征值分析,5.3.1 矩陣的基本運(yùn)算模塊,基本線性代數(shù)計(jì)算,矩陣的變換,矩陣相乘,矩陣與向量乘,IInteger SReal C Com

15、plex DDouble Z Double complex SDSingle and Double CZ Single and double complex DQDouble and Quadruple ZQDouble and quadruple complex,基本代數(shù)運(yùn)算子程序名的第一個(gè)字符的含義：,例如：SGEMM Matrix-matrix multiply, general,表示單精度的矩陣與矩陣的乘法,例5.3.1 計(jì)算一個(gè)矩陣A 與其轉(zhuǎn)置陣的乘積，這里,打開,注意B是輸出變量,編制程序：在工程中建立數(shù)據(jù)文件，在程序中打開數(shù)據(jù)文件 和結(jié)果文件，再調(diào)用矩陣乘子程序和矩陣打印程序。,

16、5.3.2 解線性方程組系統(tǒng),點(diǎn)擊線性系統(tǒng)，可看到列表,解線性系統(tǒng)的子程序名的意義：,解方法示意圖,程序名的意義：,LSARG,高精度求解,一般矩形矩陣,LFCRG,三角分解與條件數(shù),一般矩形矩陣,LSACG,矩形復(fù)數(shù)陣,高精度求解,LSVRR,實(shí)系數(shù)矩陣,奇異值分解,例5.3.2： 對系數(shù)矩陣為正定對稱陣的線性方程組,先對系數(shù)矩陣進(jìn)行Cholesky分解，再求解線性方程組 解：我們使用LCHRG子程序 功能：對對稱正定矩陣、對稱半正定矩陣進(jìn)行 Cholesky 分解。 現(xiàn)在我們來看LCHRG子程序。 Usage CALL LCHRG (N, A, LDA, PIVOT, IPVT, FAC,

17、 LDFAC),編程，先將對稱正定矩陣建立一個(gè)數(shù)據(jù)文件CHE.DAT,然后編程如下：,計(jì)算結(jié)果為：,5.3.3 求矩陣的譜分解,求矩陣的譜分解，求矩陣的特征值和特征向量。常規(guī)線性特征值系統(tǒng)求解問題是 這里為一個(gè)矩陣，而為特征值，為對應(yīng)于的一個(gè)向量。 廣義線性特征值系統(tǒng)求解問題是,常規(guī)線性特征值系統(tǒng)求解程序以 E 開頭，而廣義線性特征值系統(tǒng)求解以 GE 開頭，具體的程序分類見下列表格。,例 5.3.3 計(jì)算以下矩陣的全部特征值和特征向量,Usage CALL EVCRG (N, A, LDA, EVAL, EVEC, LDEVEC) Arguments N Order of the matrix

18、. (Input) A Floating-point array containing the matrix. (Input) LDA Leading dimension of A exactly as specified in the dimension statement of the calling program. (Input) EVAL Complex array of size N containing the eigenvalues of A in decreasing order of magnitude. (Output) EVEC Complex array containing the matrix of eigenvectors. (Output) The J-th eigenvector, correspo