回归问题——线性方程组求解的迭代法

1、第六章回归问题6.1 回归问题1.1.1 问题的引入在数理统计中，把研究对象的全体称为总体，而把组成总体的每个单元称为个体， 要了解总体的规律性，必须对其中的个体进行统计观测。但若对全部个体进行观测，这样能对总体有充分的了解，但实际上行不通，而且也不经济。所以对整体进行随机抽样观测，再根据抽样观察的结果来推断总体的性质成为一种重要的方法。许多数理统计建模的实际问题中，一个随机变量与另一个随机变量的关系不是线性关系，而是曲线关系，那么如何确定回归方程呢？下表给出了某种产品每件平均单价y（元）与批量x（件）之间的关系的一组数据，试确定y 与 x的函数关系。表 6.1.1 已知数据x20253035

2、4050606570758090y1.811.701.651.551.481.401.301.261.241.211.201.181.1.2 模型的分析先将表 6.1.1 中的数据进行曲线拟合，然后根据经过拟合的曲线形状确定回归方程的次数。用 MATLAB 做出拟合图如下，由下图知,可建立二次回归多项式模型。图 6.1.1 散点图1.1.3 模型的假设假设上表给出的数据是真实的，且以上数据是随机抽取的可以较准确地推断单位与批量的关系，假设单价与批量的函数关系是一个多项式函数，可用多项回归来建立模型。1.1.4 模型的建立根据模型的分析，可以建立多项式模型y 01x2x2,N(0, 2)，令 x

3、1 x,x2 x2，则回归方程可写成y 01x12x12,N(0, 2)，这是一个二元线性回归模型。且XTXXTY ，其中：1111111111112025303540506065707580904006259001225160025003600422549005625640081001.181.701.651.551.481.401.301.261.241.211.201.180=126.2 线性方程组迭代法概述迭代法： 即用某种极限过程逐步逼近线性方程组精确解的方法。迭代法具有需要计算机存储较少、程序设计简单、原始系数矩阵在计算过程中始终不变等优点， 但有收敛性或收敛速度的问题。迭代法是解

4、大型稀疏矩阵方程组的重要方法。迭代法能保持矩阵的稀疏性，具有计算简单、编制程序容易的优点，并在许多情况下收敛较快。故能有效地解一些高阶方程组。迭代法的基本思想是构造一串收敛到解的序列，即建立一种从已有近似解计算新的近似解的规则。由不同的计算规则得到不同的迭代法。对线性方程组AX b ， 其中， Aaij n n 非奇异矩阵，b b1 , bn T。构造其形如x Mx g 的同解方程组，其中M 为 n阶方阵，g Rn。任取初始向量x0 Rn， 代入迭代公式xk 1 Mx k g k 0,1,2, 产生向量序列 x k ，当 k 充分大时，x k 作为方程组AX b 的近似解，这就是求解线性方程组

5、的单步定长线性迭代法。M 称为迭代矩阵。6.3 迭代法6.3.1 Jacobi 迭代法a11对 Ax b，设 det(A) 0，aii0(i 1 n) (6.3.1) ，将 A改写成：0 a12a13a1n0a23a2nD L U (6.3.2)an 1,n0a21a22ana31a320an2an,n 1145 / 18又将 A分裂为：A M N ，则 (6.3.1) 等价于 Mx Nx b (6.3.3)其中 M 应选择为一个非奇异阵，并使Mz f 容易求解。对应 (6.3.3) 可构造一个迭代过程：初始向量x(0)，x(k 1) M 1Nx(k) M 1b,(k 0,1, )(6.3.4

6、)特别地，若选取M D ，则 N M A L U ，从而 (6.3.1) 化为：Dx (L U )x b可得： Jacobi 迭代公式：x(0) (初始向量)， x(k 1)Jx(k)f ，其中：J D 1(L U),f D 1b(6.3.5)J 称为 Jacobi 迭代的迭代矩阵。Jacobi 迭代的分量形式：引进记号：xkx1k , x2k ,xnk T为第 k 次近似，由(6.3.5) 有：x0x10 ,x20 ,xn0 Txik 11aiibinaijxjkj1ji, i 1 n,k 0,1,(6.3.6)Jacobi 迭代公式简单，由公式(6.3.5), (6.3.6) 可知，每迭代

7、一次只需计算一次矩阵与向量乘法，计算机中只需要两组工作单元用来保存x(k)及x(k 1)且可用x(k 1) x(k) 来控制迭代终止。由迭代计算公式可知，迭代法一个重要特征是计算过程中原来矩阵A数据始终不变。例 6.3.1 用 Jacobi 迭代法求下面线形方程组，其精确解是x* (1,2, 1,1)T：10x1 x2 2x36x1 11x2 x3 x4 252x1 x2 10x3 x4113x2x3 8x4 15解 先将 Ax b 转化为等价方程组：1x110(6 x2 2x3)x2x3x41(25 x1 x3 3x4 )11110( 11 2x1 x2 x4 )1(15 3x2 x3 )8

8、迭代公式：选取初始向量x(0) (0,0,0,0) T ，x1(k 1)x2(k 1)(k 1) x3x4(k 1)1 (6 x2(k) 2x3(k) )101 (25 x1(k)x3(k) 3x4(k) )11,k 0,1,21 ( 11 2x1(k) x(2k)x4(k)1(15 3x2(k) x3(k) )*(10)xx0.0002。经 10 次迭代解：x(10) (1.0001,1.9998, 0.9998,0.9998)T ，误差为：6.3.2 Gauss-Seidel 迭代法在 (6.3.3) 中选取 M D L(下三角阵)，则 N M A U ，从而 (6.3.1)化为等价的：(

9、D L)x Ux b(6.3.7)可得 Gauss-Seidel 迭代公式：初始向量x(0) ， xk 1 Gx k f ，其中1G DLU1(6.3.8)f DLbG 称 为 Gauss-Seidel 迭 代 矩 阵 。 G-S 迭 代 法 的 分 量 形 式 为 ： 记x kx1k ,x2k ,xnk T ，有迭代公式(6.3.9) ，Tx x1 ,x2 , xn1 i1nxi biaij xjaij xj(6.3.9)aiij 1j i 1i 1n,k 0,1,2,Gauss-Seidel 迭代法每迭代一次只需计算一次矩阵与向量的乘法，但G-S迭代法比Jacobi 迭代法有一个明显的优点

10、，那就是计算机上仅需一组工作单元用来保存x(k)分量(或x(k 1)分量)当计算出x(jk 1)就冲掉旧的分量x(jk)。由G-S迭代公式 (6.3.9) 可看出在x(k)x(k 1)的一步迭代中，计算分量xi(k 1)时利用了已计算出来的新分量x(jk 1)(j 1 i 1)。 因此， G-S迭代法可以看作是Jacobi 迭代法的一个修正。例 6.3.2 用G-S方法解下面方程组，其精确解为：x 1,2, 1,1 T，10x1 x2 2x3 6x1 11x2 x3 3x4 252x1 x2 10x3 x 4113x2 x3 8x415解 由 (6.3.9) 可得本题G-S迭代公式为：x1k

11、11 6 x2k 2x3kx2k 11 25 x1k 1 x3k 3x4k11,k 0,1,2x3k 1111 2x1k 1 x2k 1x4kx4k 11 15 3x2k 1x3k 14823经 5 次迭代得：x0(0,0,0,0) T， x51.0001, 2.0000,1.0000, 1.0000 T，x x 50.0001 。从此例可见，G-S 迭代法比Jacobi 迭代法收敛快(初始向量相同，达到同样精度，所须迭代步数较少)，但这个结论对Ax b的矩阵 A满足某些条件才是对的，甚至有这样的线性方程组，用Jacobi 方法是收敛的，而用G-S迭代法却是发散的。x1 2x2 3x3 1此例

12、用 Jacobi迭代法收敛而G-S迭代法却发如线性方程组：x1 x2 x3 12x1 2x2 x31散。 简要分析如下：1GG D L U(GG) 2 1 ，GJ0 , GJ( J1000D1L U) 。01110000特征方6.3.3 迭代法的收敛性定义 6.3.1 设有矩阵序列Ak(ai(jk)n n,(k 1,2, )及 A (aij)，如果lkim aijkaij， i, j 1 n 成立，则称Ak 收敛于 A。记为lkim AkA。1例如，矩阵序列A, A022 02Akkk10kk1 时，k 00Ak00矩阵序列收敛的概念可以用任何矩阵范数来描述。下面介绍向量、矩阵的范数定义、常用

13、的范数形式以及相关性质。定义 6.3.2 (向量范数)如果x Rn的某个实值函数N(x) |x|， (1) 正定性：|x| 0,|x| 0 x 0 (2) |cx| |c| |x|， c为实数 (3) 三角不等式：|x y| |x| |y| x.y Rn，则称N(x) |x|为 Rn上的一个向量范数(或向量。 利用三角不等式容易证明：(4) x y x y在 Rn中，三角不等式即为三角形两边长之和大于第三边长。如下图：定义633 设x=(0X2,.xn ) eRn ,定义可上四种常用向量范数如下:向量的1一范数:向量的8 范数:(3)向量的2范数:1/ f n l|x|2 = (x,x) 2

14、=|E it422XJ向量的p范数:l|x lip = 工 I Xi IIf J1p,P1,+x )易证明：它们均满足向量范数定义 6.3.1 ,且前三种向量范数均为第四种的向量范数的特例。其中(| xl lim | x|p )P例 6.3.3 设x=(-123； E 酎计算 HxILIIxILellxlbo|x|=1+2+3=6,|xb= max 1,2,3 =3解1.一|x|b =(f+22 +32 ) 2 =定理631设非负实值函数N(x)x|为上任一向量范数，则N(x)是X分量Xi X2Xn的连续函数。定理6.3.2 (范数等价性)设|卜|/以|/为川上任意两种范数，则存在常数G 6

15、>0,使得:d|x|s<|x|z<c2|x|s, VxSRn0定理6.3.3 在 山中，lim x)=x* u |x)X* IG 0 ( kT -时)其中为向量 的任一种范数。定义634 (矩阵的范数)若矩阵AE的某个非负实值函数N(A)=|A|涉足 条件：正定性：|A庐。且| A|=0J A=0(2)正齐次性:|cA|=|c| |A|, c R1(3) 三角不等式：|A B| |A| |B|(4) |AB | |A| | B |则称 N A 为上 Rn n的一个矩阵范数或模。如果将矩阵A看成Rn n向量空间中的元素，则由向量2范数定义有：n12F A | A |Fai2,

16、j。 这就是矩阵的Frobenius 范数 , 明显它满足上面定义。i,j 1在大多数与估算有关的问题中，矩阵和向量会同时参与讨论，所以希望引进一种矩阵的范数。它和向量范数相联系且和向量范数是相容的，即：|Ax| |A| |x|, x Rn, A Rn n定义 6.3.5 (矩阵的算子范数)设x Rn, A Rn n给出一种向量范数| x| (如1,2或) ，相应地定义一种矩阵的非负函数：|A| max |Ax| (最大|x| 0 |x|比值) 易验证： | A| 满足前一矩阵范数的定义。因此 | A| 是 Rn n上的一个范数，称之为 A的算子范数。定理 6.3.4 设 |x| 是 Rn上的

17、一个向量范数，则|A| 是 Rn n上的矩阵范数，且满足相容条件：|Ax| |A| |x| 。定理 6.3.5 设 x Rn,A Rn n则：n(1) |A| max |aij | (称为 A的行范数) 1in j1 n(2) | A |1 max| aij | (称为 A的列范数)1jni1(3) |A|2max ATA(称为 A的2范数)其中 max AT A 表示AT A 的最大特征值。在计算上，(1) ， (2) 比较容易，而(3) 比较困难，但(3) 在理论分析上十分有用。定理 6.3.6lkim AkA充要条件是Ak A 0, k ，我们要考虑迭代法收敛 性 ， 即 要 研 究 B

18、在 什 么 条 件 下 使 误 差 向 量 趋 于 零 向 量 。 即k x k xBk 00, k定理 6.3.7 设 B (bij )n n， 则 Bk0(零矩阵)， (k )的充要条件是：(B) 1 。定理 6.3.8 (迭代法基本定理)设有方程组x Bx f ， 对于任意初始向量x(0)及任意 f ， 解此方程组的迭代法(即x(k 1)Bx(k) f ) 收敛的充要条件是：(B) 1 。此定理应用时要计算谱半径，不太方便。将其修改成如下的形式：定理 6.3.9 (迭代法收敛的充分条件)设有方程组x Bx f ，且 x(k) 为由迭代法产生的序列x(k 1)Bx(k)f ， x(0)为初

19、始任意向量。若迭代矩阵B 的某种范 数满足 B q 1 ，则：(1) x(k) 收敛于方程组I B x f 的唯一解x* 。(2) x*x(k) 11xk( 1) xk。()1qk(3) 误差估计：x x k q x1 x 0 。1q证明 (1) 由于 B q 1 ，即 I B可逆，有唯一解。设为x ， x Bx f 。引进误 差 向 量 ek xk x ， 即 得 误 差 ek 的 递 推 公 式 ： ek 1 BekBk1e0 ,k 0,1, ， 于是 ek Bke0qke00，k时 e0x0 x。(2) 由迭 代 公 式 有 xk 1 xk B xk xk 1 ， 又 ek1 B ek

20、则xk 1 xk x xk x xk1 x xk x xk 11 q x xk ，即10x1 x2 2x36例 6.3.4 考察用 Jacobi 方法求解线性方程组x1 11x2 x3 x4 25 的收敛2x1 x2 10x3 x4113x2x3 8x4 15性。10121 111解A21100310318101110012800100103102013010则 Jacobi 迭代矩阵为：J D1L U111210011021000111311110011038180J max 3 , 5 , 4 , 41 1. 故解此方程组的Jacobi 方法收敛。10111082k xk xBk 0 ，

21、设 B 有 n 个线性无关的特征n向量u1, u2un，相应的特征值为1，2 n，由 kaiui (线性表示，基)i1nnkBk 0ai Bkuiai ikui 。 可 以 看 出 当 (B) 1 愈 小 时 ，i1i1ik 0 i( 1n 。 (k )愈快，即(k)0愈快，故可用(B)来刻划迭代法的收敛快慢。现在来确定迭代次数k ，使 (B) k 10 s，取对数得：k sln10 。ln (B)定义 6.3.6 称 R(B) ln ( B)为迭代法的收敛速度。由此可见，(B) 1 愈小，ln (B) 愈大，则所需的迭代次数就越少。而 n 较大时，矩阵特征值计算比较困难，基本定理的条件比较难

22、验证，所以最好建立与矩阵元素直接有关的条件来判断迭代法的收敛性。由于(B)vB v， 所以可用B v作为(B)v上界的一种估计，这样的结果即为迭代法收敛的充分性条件(如上面的定理)。由这个充分性条件的结果(3) 可见， B q 1 越小，收敛越快。另外，由其结果(2) 可知：若x(k)x(k 1) 则x* x(k) q 0。因此可以用1qxk xk1 作为控制迭代终止的条件。但要注意当q 1 时， q 较大，尽管x(k) x(k 1) 很小，却有可能误差向量模q1(1) x*x(k) 很大，迭代法收敛很慢。而且此充分性条件的结果(3) 可以用来事先确定需要迭代多少次才能保证(k)。定理 6.3

23、.10 解方程组Ax b的 G-S迭代法收敛的充分必要条件是(G) 1 ， G为 G-S迭代矩阵。定义 6.3.7 (对角占优矩阵)设A (aij)n nn(2) 若aiiaij0(i 1 n) ，即 A的每一行对角元素绝对值严格大于同行其j1 ji他元素绝对值之和。则称A为严格对角占优矩阵。n(3) 若aiiaij 0(i 1 n) ，且至少有一个不等式严格成立，则称A为弱对j1 ji角占优矩阵。定义 6.3.8 (可约与不可约矩阵)设A (aij)n n，当 n 2时，如果存在n阶排列AA矩阵 p使pTAp1112 成立。其中A11 为r 阶子矩阵，A22为n r 阶子矩阵0A22( 1

24、r n) ， 则称矩阵A为可约的。若不存在排列矩阵p 使上式成立，则称 A为不可约矩阵。当 A为可约矩阵，则 Ax b可经过若干行列重排化为两个低阶方程yd组求解。 事实上， 由 Ax b化为：PTAP(PTx) PTb， 设 y PTxy1 ,PTbd1 ，y2d2其中y1,d1为 r 维向量。于是，求解Ax b化为求解：A11y1 A12y2 d1 ，另外 A为可约矩阵的充要A22y2 d2条件：存在指标集J 1,2, ,n ,J 使akj 0， k J， j J.定理 6.3.11 (对角占优定理)：若 A (aij )n n为严格对角占优阵或为不可约弱对角占优阵，则A 是非奇异阵。证明

25、 (1) A为严格对角占优阵，采用反证法。若 detA 0， 则 Ax 0有非零的解，设 为 x (x1x2xn )T0 。设xkm axxi1in ixi 0齐 次 方 程 中 的 第 k个 方 程 得 ：nakix0，则可得：ki ji1akk xknakjxjj1jknakjj1jkxjnxkj1 jk(6.3.10)n即有：akkakj 这与严格对角占优阵矛盾，故det A 0。j1 jk(2) A为不可约弱对角占优阵，采用反证法。n设 x 0, x (x1xn)T 使Ax 0。并令 m 使ammamj .(6.3.11) (由弱对角占j1 jm优定义)成立。再定义下标集合：J k x

26、kxi , i 1 n, xkxj .对某个j在 (6.3.10) 中取 k m，将导致与(6.3.11) 得矛盾。故J (空集合)。对任意nk J ， 有akjka j xk/ x (由 (6.3.10) ) ， 由此可见，当xkxj 时有 akj 0。j1jk否则，上式就与A为弱对角占优阵矛盾。但对任意k J 和 j J ，必有 xk xj ，因而akj0,k J,j J 从而 A为可约矩阵，这与A 为不可约矩阵假设矛盾。定理 6.3.12 若 A Rn n为严格对角占优矩阵或为不可约对角占优矩阵，则对任意的初始向量x(0) ，方程组Ax b的Jacobi 迭代法和G-S迭代法均收敛。且G

27、-S迭代法比Jacobi 迭代法收敛速度快。证明 设 A为不可约对角占优矩阵，先证明G-S迭代法收敛。由弱对角占优阵假设知aii 0,(i 1 n)，而 G-S 迭代矩阵为G (D L) 1U ，又 det(D L) 1 0， det( I G) det( I (D L) 1U) det(D L) 1， det( (D L) U) 0即 det( (d l) U ) 0.，记 C (D L) U 。则：a11a12a13a1nca21a22a23a2nan1an2an3ann下面来证明当1 时， detC 0， 这是因为A为不可约阵，则 C也不可约，由 A为弱对角矩阵，可得：i1n(当1)ci

28、iaiiaijaijcijj1ji1ji且至少有一个不可约等式严格成立，这表明当1 时， C 为不可约弱对角占优矩阵，于是由前一个定理可知当1 时， detC 0，这一结论表明detC0的根 满足：1 ，即 (G) 1 ，故G-S法收敛。在同一条件下，对于Jacobi 方法( J D 1(L U )完全类似于上可证。当A为严格对角占优阵时，证明方法完全类似。6.3.4 超松弛代法逐次超松弛迭代法(Successive Over Relaxation Method. 简称为SOR法)是 Gauss-Seidel 法的一种加速方法。这是解大型矩阵方程组的有效方法之一，具有计算公式简单、程序设计容易

29、、占用计算机内存较少等优点，但需要选择好的加速因子，即最佳的加速因子。对 Ax b(6.3.12) 其中 A Rn n为奇异矩阵，且设aii0.(i 1 n)，分解：A D L U(6.3.13)设 已知第 k 次 迭 代 向 量 x(k) ， 及 第 k 1 次 迭 代 向 量 xi(k1 的)分 量x(jk1() j1 , 2 ,i ，现在来计算分量： 1 )先用 G-S迭代法求出辅助量xi(k 1)(预测)i1nxi(k 1)(biaijx(jk 1)aijx(jk),i 1 n (6.3.14)aiij 1j i 1再取xi(k 1)为 xi(k)与 xi(k 1)的某种平均值(即加权

30、平均)，即xi(k 1)(1)xi(k)xi(k1)xi(k)(xi(k1)xi(k)(6.3.15) (校正)将 (6.3.14) 代入 (6.3.15) 即得解 Ax b的逐次超松弛迭代公式：162 / 18(k1)(k)i1 (k1) n (k)xixi(biaij xjaijxj )aiij 1j ix(k)(x1(k),x2(k),xn(k),(k 0.1, , i 1 n)其中称 为松弛因子，或写为：xi(k 1)xi(k)xii1n(k 1)(k)xi(biaij xjaij xj )aiij 1j i(k 0,1, , i 1,2,n)(6.3.16)(6.3.17)显然 1

31、时，解 (6.3.12) 的 SOR法即为Gauss-Seidel 迭代法。SOR 法中每迭代一次，主要计算量是计算一次矩阵与向量乘法。由(6.3.16)可见在计算机上用SOR法解方程组只需一组工作单位元，以便存放近似解。迭代计算时，可用p max ximax xi(k 1) xi(k)来控制迭代终止。当 1 时，称 (6.3.16) 为低松弛法；当1 时，称 (6.3.16) 为超松弛法。例 6.3.5 用 SOR法解：14111141其精确解为x*( 1, 1, 1, 1)T解 取 x(0)0 ，则SOR迭代公式为：x1(k 1)x1(k)(14x1(k)x2(k)x3(k)x4(k)/4

32、x2(k1)x2(k)(1x1(k1)4x2(k)x3(k)x4(k)/4x3(k1)x3(k)(1x1(k1)x2(k1)4x3(k)x4(k)/4x4(k1)x4(k)(1x1(k1)x2(k1)x3(k1)4x4(k)/4取 1.3 ，第 11 次迭代结果为：x(11) ( 0.99999646, 1.00000310, 0.99999953, 0.99999912)T(11)x(11)x*0.46 10 5松弛因子满足误差x k x* 10 5 的迭代次数1.0221.117对 取其它值，迭代次数如下表，由此例可见，松弛因子选择得好，会使SOR迭代的收敛大大加速。本例中，1.3是最佳松

33、弛因子。表 6.3.1 结果表1.2121.3*11* 最少迭代次数1.4141.5171.6231.7331.8531.9109下面考察SOR迭代公式的矩阵形式，由(6.3.16) 可改写为：i1n(k 1)(k)(k 1)(k)aii xi(1)aii xi(biaij xjaij xj ), i 1 n (6.3.17)j1ji1由 A D L U ，则：Dx( k 1) (b Lx(k 1) Ux(k) (1 )Dx(k)即D L)x(k1) (1)DU )x(k)b，由于对任意，(D L)均为奇异阵(设aii0.i 1 n, 而 L 为下三角阵，且对角线元素为0)则：x(k 1)(D

34、L) 1 (1)DU x(k) (DL)1b.因此，若aii0,i 1 n ，则SOR迭代公式为：x(k 1)L x(k)f(6.3.18)其中 L (D L) 1 (1 )D U , f (D L) 1b， 称 L 为SOR方法的迭代矩阵，应用关于迭代法的收敛性定理于(6.3.18) 可得：定理 6.3.13 对 Ax b，且aii 0(i 1 n) 则解方程组的充要条件是：(L ) 1。引进超松弛法的想法是希望能选择松弛因子使迭代过程(6.3.16) 收敛较快，也就是应选择因子使 ( L *)min ( L *) 。下面考虑对于方程组(6.3.12) ( aii 0,i 1 n) ，超松弛因子在什么范围内取值才可能收敛。定理 6.3.14 (必要条件)对Ax b， ( aii,0 i1 n)的SOR方法若收敛，则：02。证明 由 SOR法收敛及上定理，知( L )1 ，设 L 的特征值为1，2， ， n则：1det(L )1 2 n ( (L )n，即 det(L )n (L ) 1，而det( L ) det( D L) 1)det(1 )D U) (1)n ， 因 此 11 ， 即 ：02。此结果表明要使SOR法收敛，松弛因子必须在 0,2 中。那么，反过来，若选取 在 0,2 中，SOR法是否一定收敛呢？定理 6.3.15 (