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

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)来刻划迭代法的收敛快慢。现在来

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)