1、第六章回歸問題線性方程組求解的迭代法6.1回歸問題6.1.1問題的引入在數理統計中，把研究對象的全體稱為總體，而把組成總體的每個單元稱為個體，要了解總體的規律性，必須對其中的個體進行統計觀測。但若對全部個體進行觀測，這樣能對總體有充分的了解，但實際上行不通，而且也不經濟。所以對整體進行隨機抽樣觀測， 再根據抽樣觀察的結果來推斷總體的性質成為一種重要的方法。 許多數理統計建模的實際問題中， 一個隨機變量與另一個隨機變量的關系不是線性關系， 而是曲線關系，那么如何確定回歸方程呢？下表給出了某種產品每件平均單價 y（元）與批量 x（件）之間的關系的一組數據，試確定 y 與 x 的函數關系表6.1.1

2、已知數據x202530354050606570758090y1.811.701.651.551.481.401.301.261.241.211.201.186.1.2模型的分析先將表 6.1.1 中的數據進行曲線擬合，然后根據經過擬合的曲線形狀確定回歸方程的次數。用 MATLAB 做出擬合圖如下，由下圖知，可建立二次回歸多項式模型。圖6.1.1散點圖6.1.3模型的假設假設上表給出的數據是真實的， 且以上數據是隨機抽取的可以較準確地推斷單位與批量的關系， 假設單價與批量的函數關系是一個多項式函數， 可用多項回歸來建立模型。6.1.4模型的建立根據模型的分析，可以建立多項式模型 y=瓦+Px+%

3、x2+小 N(0,62),令 x1=x,x2=x2,則回歸方程可寫成 y=P0+P1X1+P2x2+%1N(0,62),這是6.2線性方程組迭代法概述迭代法：即用某種極限過程逐步逼近線性方程組精確解的方法。迭代法具有需要計算機存儲較少、程序設計簡單、原始系數矩陣在計算過程中始終不變等優點，但有收斂性或收斂速度的問題。 迭代法是解大型稀疏矩陣方程組的重要方法。 迭代法能保持矩陣的稀疏性，具有計算簡單、編制程序容易的優點，并在許多情況下收斂較快。故能有效地解一些高階方程組。迭代法的基本思想是構造一審收斂到解的序列，即建立一種從已有近似解計算新的近似解的規則。由不同的計算規則得到不同的迭代法。對線性

4、方程組 AX=b,其中，A=(aj幾非奇異矩陣，b=(b1,III,bnfo構造其形如 x=Mx+g 的同解方程組，其中 M 為 n 階方陣，gWRn。任取初始向量 x/Rn,代入迭代公式 x(k4t)=Mx2+g(k=0,1,2j|)產生向量序列x(k),當 k 充分大時，x(k 府為方程組 AX=b 的近似解，這就是求解線性方程組的單步定長線性迭代法。M 稱為迭代矩陣。個二元線性回歸模型。(XTX)P=XTY,其中:一11111！111112025303540506065707580904006259001225160025003600422549005625640081001.181.7

5、01.651.551.481.401.301.261.241.211.201.1816.3迭代法6.3.1 Jacobi迭代法又將 A 分裂為：A=MN,則(6.3.1)等價于 Mx=Nx+b(6.3.3)其中 M 應選擇為一個非奇異陣，并使 Mz=f 容易求解。對應(6.3.3)可構造一個迭代過程：初始向量 x(0),x(k1)-MNx(k)M,b,(k-0,1,HI)(6.3.4)特別地，若選取 M=D,則 N=M-A=L+U,從而(6.3.1)化為：Dx=(LU)xb可得：Jacobi 迭代公式：x(0)(初始向量)，x(k=Jx(k)+f,其中：J=D(LU),f=D,b(6.3.5)

6、J 稱為 Jacobi 迭代的迭代矩陣。Jacobi 迭代的分量形式：引進記號：xN=(x!k)x2kl川婿日為第 k 次近似，由(6.3.5)有:x0=XixHIxn0TJacobi 迭代公式簡單，由公式(6.3.5),(6.3.6)可知，每迭代一次只需計算一次矩陣與向量乘法，計算機中只需要兩組工作單元用來保存 x(k)及 x(k+)且可用|x(kWx(k)|資8 來控制迭代終止。由迭代計算公式可知，迭代法一個重要特征是計算過程中原來矩陣 A 數據始終不變。例 6.3.1 用 Jacobi 迭代法求下面線形方程組，其精確解是 x*=(1,2,-1,1)T：a11對 Ax=b，設 det(A)

7、#0,a.=0(i=11n)(6.3.1)，將 A 改寫成：I0-a12-al3|-aln|0-323|-32n一,kD-L-U(632)1aiJn-0J-321a?2-a31-as20aIN-ai.nlxik1=-abin-工aM)j1j-i：,i=1n,k=0,1,|(6.3.6)10X1X2+2&=6-K+11x2-X3+X4=252KX2+10&X4=113X2X3+8X415解先將 Ax=b 轉化為等價方程組:1sX1=10(6X2-2X3)迭代公式：選取初始向量 x(0)=(0,0,0,0)T,經 10 次迭代解：x(10)=(1.0001,1.9998,-0.9998,0.999

8、8)T,誤差為：|6.3.2 Gauss-Seidel迭代法在(6.3.3)中選取 M=D-L(下三角陣)，則 N=M-A=U,從而(6.3.1)化為等價的：可得 Gauss-Seidel 迭代公式：初始向量 x(0),x(k+)=Gx(k)+f,其中1G-ID-LU1LiD-LbG 稱為 Gauss-Seidel 迭代矩陣。G-S 迭代法的分量形式為：記X2X3X41(25x1x3-3x4)111，一 c、-10(-11-2x1x2X4)1.c、X1(k1)(k1)x3x4k1)(6+x2(k)-2x3(k)10(25+x，+x3k)-3x4k)?,k=0,1,2|(112x1(k)+x2k

9、)+x4k)=1=(153x2k)+x3k)*(10)x-x/0.0002。(D-L)x=Uxb(6.3.7)(6.3.8)乂 6=口，)x$)用 x，)T,有迭代公式(6.3.9),T-TTxf)=(寸)?孰|X)1inxf*)=bi-ZaM*)-a/，)(639)aiij4I 書ji=1n,k=0,12 川Gauss-Seidel 迭代法每迭代一次只需計算一次矩陣與向量的乘法，但 G-S迭代法比 Jacobi 迭代法有一個明顯的優點，那就是計算機上僅需一組工作單元用來保存 x(k)分量(或 x(k*分量)當計算出 xjk用就沖掉舊的分量 x(k)。由 G-S 迭代公式(6.3.9)可看出在

10、 x(k)Tx(k41)的一步迭代中，計算分量 x，*)時利用了已計算出來的新分量 xjk+)(j=1|_i-1)。因此，G-S 迭代法可以看作是 Jacobi 迭代法的一個修正。例 6.3.2 用 G-S 方法解下面方程組，其精確解為：x*=(1,2,1,1f，10 x1-x2+2x3=6_為+11x2x3+3x4=2512x1-x210 x3-x4=-113x2-x38x4=15解由(6.3.9)可得本題 G-S 迭代公式為：eq2kkbx2k*)=工(25+xM)+x3k)-3x4khJ,k=0,1,2111x)=.11_2x1(k*)+x2k”x4k)x4ks=1(153xJk+)+/

11、*)L8經 5 次迭代得：x()=(0,0,0,0)T,x(5)=(1.0001,2.0000,-1.0000,1.0000T,|x*-xf 匕生 0.0001。從此例可見，G-S 迭代法比 Jacobi 迭代法收斂快(初始向量相同，達到同樣精度，所須迭代步數較少)，但這個結論對 Ax=b 的矩陣 A 滿足某些條件才是對的，甚至有這樣的線性方程組，用 Jacobi 方法是收斂的，而用 G-S 迭代法卻是發散的。x12x2-3x3=1如線性方程組：x1x2x3=1J2x12x2x3=1散。簡要分析如下:矩陣序列收斂的概念可以用任何矩陣范數來描述。下面介紹向量、矩陣的范數定義、常用的范數形式以及相

12、關性質。定義 6.3.2(向量范數)如果 xWRn的某個實值函數 N(x)引 x|,(1)正定性:|x|盧0,|x|=0Hx=0(2)11cxi 付 c|,|x|,c 為實數(3)三角不等式:|x+y 舊|x|+|y|x.ywRn,則稱 N(x)三|x|為 Rn上的一個向量范數(或向量的模)。利用三角不等式容易證明:在 Rn中，三角不等式即為三角形兩邊長之和大于第三邊長。如下圖:,止匕例用 Jacobi迭代法收斂而 G-S 迭代法卻發jGGh：D-LU00-2003-101-1Y001-20Y00人0-2003-10001,GJ-1-2-20-23-10特征方=0：；Gj(=J6.3.3迭代法

13、的收斂性定義 6.3.1 設有矩陣序列入=(a；k):(k=1,2,|)及 A=(a3如果kmaf)=純,(i,j=1n)成立，則稱入收斂于A。記為叩 A，一，一一1例如，矩陣序列 A=|,A：022、=_0叫 l!AkkkJ-10kk,當0,使得:d|x|s|x|z的一個矩陣范數或模。如果將矩陣 A 看成 Rnn向量空間中的元素，則由向量 2范數定義有：/n丫2F（A）|A|F=工 aj。這就是矩陣的 Frobenius 范數，明顯它滿足上面定義。J4）在大多數與估算有關的問題中，矩陣和向量會同時參與討論，所以希望引進一種矩陣的范數。它和向量范數相聯系且和向量范數是相容的，即：11Ax|_|

14、A|x|,-xRn,-ARnn定義 6.3.5（矩陣的算子范數）設 xwRn,A 三 Rn沏給出一種向量范數|x|1V（如（如 v=1,2 或+8），相應地定義一種矩陣的非負函數：11A|=max11Axi1（最大.網-0|x|L比值）易驗證：|A1 滿足前一矩陣范數的定義。因此|A11V是 Rn刈上的一個范數，稱之為 A 的算子范數。定理 6.3.4 設|x|v 是 Rn上的一個向量范數，則|A|是 Rn沏上的矩陣范數，且滿足相容條件:|Ax|A|J|x|L。定理 6.3.5 設 xwRn,Aw 田坨則：n|Ak=max|4|（稱為 A 的行范數）1:inj1n|A|1=maxZaj|（稱為

15、 A 的列范數）11A|2=Gmax（ATA）（稱為 A 的 2范數）其中九max（ATA 宸示（ATA 炳最大特征值。在計算上，（1）,（2）比較容易，而（3）比較困難，但在理論分析上十分有用。定理 6.3.6pmA=A 充要條件是 IIA.-AH0,（kT8）,我們要考慮迭代法收斂性，即要研究B 在什么條件下使誤差向量趨于零向量。即/)=x(k)-X*=Bk40)T0,(kT0定理 6.3.7 設 B=（bj）nM,則 B；0（零矩陣），（kTg）的充要條件是：P（B）1。定理 6.3.8（迭代法基本定理）設有方程組 x=Bx+f,對于任意初始向量 x（0）及任意 f，解此方程組的迭代法（

16、即 x（k）=Bx 的+f）收斂的充要條件是：P（B）10此定理應用時要計算譜半徑，不太方便。將其修改成如下的形式：定理 6.3.9（迭代法收斂的充分條件）設有方程組 x=Bx+f,且x（k）為由迭代法產生的序列 x（E）=Bx（k）+f,x為初始任意向量。若迭代矩陣 B 的某種范數滿足怛卜 q1,則：k（k）收斂于方程組（IB）x=f 的唯一解 x*。1*廣曉1-qk誤差估計：卜沖.xC 卜網 LxH。1-q證明（1）由于同=q|x*-xCj-|x*-x（k|（1-qj|Z-x（kj,即10 x1x2+2x3=6例 6.3.4 考察用 Jacobi 方法求解線性方程組+11%-=25的收斂2

17、x1-x2+10 x3-M=T13x2-&8x4=15利用此遞推式即可得誤差估計則 Jacobi 迭代矩陣為:%0_%0010%-%1%00%0一%X0.Jll=maxI,41=11.故解此方程組的 Jacobi 方法收斂。10111082卜面考察迭代法的收斂速度。sCUxC)-x*=Bksf),設 B 有 n 個線性無關的特征n向量312川/，相應的特征值為%,即|缸,由虱k)=aiUi(線性表示,基)ynn8C)=Bk82=aBkui=Zaiui。可以看出當 P(B)1 愈小時,i1i1%kT0i(=1n(kT 叼愈快，即 w(k)T0 愈快，故可用 P(B)來刻劃迭代法的收斂快慢?，F在來

18、確定迭代次數 k,使 IP(B)手10,取對數得：k 之sln10。-ln7(B)定義 6.3.6 稱 R(B)=-lnP(B)為迭代法的收斂速度。由此可見，P(B)1 愈小,-lnP(B)愈大，則所需的迭代次數就越少。而 n 較大時，矩陣特征值計算比較困難，基本定理的條件比較難驗證，所以最好建立與矩陣元素直接有關的條件來判斷迭代法的收斂性。由于 P(B)v|B|v,所以可用|B|v作為 P(B)v上界的一種估計，這樣的結果即為迭代法收斂的充分性條件(如上面的定理)。由這個充分性條件的結果(3)可見，|B|=q1 越小，收斂越快。另外，由其結果為控制迭代終止的條件。性。10-12刎-111-1

19、解 A=2-110-03-1013-18101110一 01-1-28一一01-010010-I-310-201-3010111-%00可知：若卜-x(k,)|8 則|x*-x(k)|=。因此可以用1-qx(k)_x(k“w 作但要注意當 q時，旦較大，盡管 x(k)-x(kT|很小，卻有可能誤差向量模 q-111|(k)P|x*-x(k)|很大，迭代法收斂很慢。而且此充分性條件的結果(3)可以用來事先確定需要迭代多少次才能保證|8(k)|&o定理 6.3.10 解方程組 Ax=b 的 G-S 迭代法收斂的充分必要條件是 P(G)0(i=1n),即 A 的每一行對角元素絕對值嚴格大于同行其jm

20、j小他元素絕對值之和。則稱 A 為嚴格對角占優矩陣。n(2)若為 a0(i=1n),且至少有一個不等式嚴格成立，則稱 A 為弱對j1j角占優矩陣。定義 6.3.8(可約與不可約矩陣)設 A=(aij)nM,當 n*2 時，如果存在 n 階排列A,A。矩陣 p 使 pTAp=七成立。其中A1為 r 階子矩陣，A22為 n-r 階子矩陣 J3(1Wrn),則稱矩陣 A 為可約的。若不存在排列矩陣 p 使上式成立，則稱 A 為不可約矩陣。當 A 為可約矩陣，則 Ax=b 可經過若干行列重排化為兩個低階方程fdC組求解。事實上，由 Ax=b 化為：PTAP(PTx)=PTb,設 y=PTx=y1,PT

21、b=d1,Zamj.(6.3.11)(由弱對角占j4j-m優定義)成立。再定義下標集合：J=k|xk|耳為，i=1|_n,xk,|xj.對某個 j在(6.3.10)中取 k=m,將導致與(6.3.11)得矛盾。故 J0 巾(空集合)。對任意nkwJ,有 ak|jkajxk/x(由(6.3.10),由此可見，當xjaxj時有 akj=0。j1j水否則，上式就與 A 為弱對角占優陣矛盾。但對任意 kwJ 和 j 皂 J,必有 xk|xj,因而 akj=0,kwJ,j 正 J 從而 A 為可約矩陣，這與 A 為不可約矩陣假設矛盾。定理 6.3.12 若 AWR.為嚴格對角占優矩陣或為不可約對角占優矩

22、陣，則對任意的初始向量 x(0),方程組 Ax=b 的 Jacobi 迭代法和 G-S 迭代法均收斂。且 G-S 迭代法比Jacobi 迭代法收斂速度快。證明設 A 為不可約對角占優矩陣，先證明 G-S 迭代法收斂。由弱對角占優陣假設知冉#0,(i=1n),而 G-S 迭代矩陣為 G=(DL),U,又 det(DL)#0,det(九IG)=det(九 I_(DL),U)=det(DL)，,det(九(DL)U)=0 即det(Md-l)-U)=0.,記 C=“DL)U。貝 U:a11a12a13a1nc=Ka21+7a22a23ia2n44仁兒an1九an2aan3一兒ann；下面來證明當惘之

23、 1 時，detC#0,這是因為 A 為不可約陣，則 C 也不可約，由 A 為弱對角矩陣，可得：i1n設為 x=(X1X2I11Xn)T0。設 Xk=maxxi1i：nxTO由齊次方程中的第 k 個方程得:akkxkn、ajjij水n-2akjj1j，kxjnxkSj1j：k(6.3.10)(當囚之 1)Cii=囚同|空，同+Z 同=Cijj1j工ij在且至少有一個不可約等式嚴格成立，這表明當出之 1 時，C 為不可約弱對角占優矩陣，于是由前一個定理可知當同上 1 時，detC#0,這一結論表明 detC=0 的根兒滿足：困1,即 P(G)1,故 G-S 法收斂。在同一條件下，對于 Jacob

24、i 方法(J=D(L+U)完全類似于上可證。當 A 為嚴格對角占優陣時，證明方法完全類似。6.3.4超松弛代法逐 次 超 松 弛 迭 代 法 (SuccessiveOverRelaxationMethod. 簡 稱 為 SORfe) 是Gauss-Seidel 法的一種加速方法。這是解大型矩陣方程組的有效方法之一，具有計算公式簡單、程序設計容易、占用計算機內存較少等優點，但需要選擇好的加速因子，即最佳的加速因子。對 Ax=b(6.3.12)其中 AwRn將為奇異矩陣，且設 a.#0.(i=1 用 n),分解：A=D-L-U(6.3.13)設已知第 k 次迭代向量 x(k),及第 k+1 次迭代

25、向量,-1的分量01,2,4，郵在來計算分量:先用 G-S 迭代法求出輔助量小池(預測)i1nWi(k(bi-兩 x(k1)-aijx(k),i=1n(6.3.14)aiij1jd1再取 x，*)為 x(k)與#(k41)的某種平均值(即加權平均)，即x(k41)=(1-o)x(k)+xfk+)=x(k)+切(幫一 x(k)(6.3.15)(校正)將(6.3.14)代入(6.3.15)即得解 Ax=b 的逐次超松弛迭代公式：(k1)(k),小;(k1);(k)、x=Xi(biJaijXj-、aijXj)aiij1jix(k)=(姆,x2k),染),(k=0.1,， i=1n),其中稱co 為松

26、弛因子，或寫為：乂卜的二乂，)十%iAnCO(k4)_(k)、X=(h-ZajXj-ZajXj)aiij二j(k=0,1，i=1,2,n)(6.3.16)(6.3.17)顯然 6=1 時，解(6.3.12)的 SOFfe 即為 Gauss-Seidel 迭代法。SOR 法中每迭代一次，主要計算量是計算一次矩陣與向量乘法。由(6.3.16)可見在計算機上用 SOFfe 解方程組只需一組工作單位元，以便存放近似解。迭代計算時，可用 pd=maX|AX=maXX，*，X(k)8 來控制迭代終止當缶1 時，稱(6.3.16)為低松弛法；當缶下 1 時，稱(6.3.16)為超松弛法例 6.3.5 用 S

27、ORt 解：1-41111-41其精確解為 X=(-1,-1,-1,-1),解取 X(0)=0,則 SO 砂代公式為:.娟書)=X1缶(1+4X1(k)-X2k)-X3k)-才)/4jX=X2k)/(LX 尸)+4X2k)X3k)X4k)/4X3k*=X3k)f(1-婿地尸)十 4X3k)-X4k)/4、x 尸)=x4k)(1x(k*)_xk)_x”)+4x4k)/4取 0=1.3,第 11 次迭代結果為：x(11)=(-0.99999646,-1.00000310,-0.99999953,-0.99999912(11)|2=|x(11)-x*|20.4610松弛因子滿足誤差 x(k)-xW1

28、03 的迭代次數1.0221.117對缶取其它值，迭代次數如下表，由此例可見，松弛因子選擇得好，會使SO 砒代的收斂大大加速。本例中，=1.3 是最佳松弛因子。表6.3.1結果表1.2121.3*11*最少迭代次數1.4141.5171.6231.7331.8531.9109下面考察 SOR4 代公式的矩陣形式，由(6.3.16)可改寫為：i1n(k1)(k)(k1)(k八aiiXi=(1 一)aiiX-(b-ajXj-、a。7),i=1n(6.3.17)j1j41由人=0LU,貝 u：Dx(k%=0(b+Lx(kd1)+Ux(k)+(1co)Dx(k)即(D-coL)x(k+)=(1-co)

29、D+coU)x(k)+ob,由于對任意切,(D。L)均為奇異陣(設 a.=0.i=1n,而 L 為下三角陣，且對角線元素為 0)則：x(k1)-(D-L)J1(1-)DUlx(k)(D-L)Jb.因此，若 aii=0,i=1n,則 SORI 代公式為：x(k1)=Lx 的 f(6.3.18)其中 L3=(DcoL)，1(1co)D 十切 U，f=co(DcoL)b,稱 q 為 SO 時法的迭代矩陣，應用關于迭代法的收斂性定理于(6.3.18)可得：定理 6.3.13 又 tAx=b,且%=0(i=1r)則解方程組的充要條件是：P(LJ10引進超松弛法的想法是希望能選擇松弛因子使迭代過程(6.3

30、.16)收斂較快，也就是應選擇因子缶使P(L*)=minP(L*)。000下面考慮對于方程組(6.3.12)(aii0,i=1n),超松弛因子在什么范圍內取值才可能收斂。定理 6.3.14(必要條件)對 Ax=b,(a.#0n)的 SO 時法若收斂，則：02。證明由 SORa 收斂及上定理，知 P(LQ1,設 q 的特征值為%，4 則：1det(L4|=九%以以心。)即det(LrP(/)1,而det(L=det(DcoL)-)det(1-)D+coU)=(1-)n,因止匕 1切|1,即：02。此結果表明要使 SORfe 收斂，松弛因子與必須在(0,2)中。那么，反過來，若選取缶在(0,2)中，SORt 是否一定收斂呢？定理 6.3.15(充分條件)若 A 為對稱正定矩陣， 且 0與2,則解(6.3.12)的 SORfe 必收斂。證明若能證明在定理條件下，對 L.、的任一特征值入有:-1兒1,則定理得證。事實上，設 y 為對應九的 L,、的特征向量。即：WLj=Ky,y=(y1y2，yn)T#0,即(D8L),(1-co)D+coU)y=ay。亦即：(1-o)D+8U)y=MD-6L)y,為找九的表達式，考慮內積：(10)