科学计算与数学建模第六章.pptx

第六章 回归问题 线性方程组求解的迭代法回归问题回归问题6.1线性方程组迭代法概述线性方程组迭代法概述6.2迭代法迭代法6.3关于回归模型的求解关于回归模型的求解6.41点击进入点击进入点击进入点击进入第1页/共93页 6.1 回归问题2点击进入点击进入点击进入点击进入第2页/共93页6.1.1 问题的导入 在数理统计中，把研究对象的全体称为总体，而把组成总体的每个单元称为个体，要了解总体的规律性，必须对其中的个体进行统计观测。但若对全部个体进行观测，这样能对总体有充分的了解，但实际上行不通，而且也不经济。所以对整体进行随机抽样观测，再根据抽样观察的结果来推断总体的性质成为一种重要的方法。许多数理统计建模的实际问题中，一个随机变量与另一个随机变量的关系不是线性关系，而是曲线关系，那么如何确定回归方程呢？3第3页/共93页x202530354050606570758090y1.811.701.651.551.481.401.301.261.241.211.201.18下表给出了某中产品每件平均单价下表给出了某中产品每件平均单价 y y（元）与批量（元）与批量 x x（件）之间的（件）之间的关系的一组数据，试确定关系的一组数据，试确定 y y与与 x x的函数关系的函数关系。表表6.1.1 6.1.1 已知数据已知数据4第4页/共93页 先将表6.1.16.1.1中的数据进行曲线拟合，然后经过拟合的曲线形状确定回归方程的次数。用MATLABMATLAB做出拟合图如下，由下图知,可建立二次回归多项式模型。5第5页/共93页6图图6.1.1 6.1.1 散点图散点图第6页/共93页6.1.3 模型的假设 假设上表给出的数据是真实的，且以上数据是随机抽取的可以较准确的推断单位与批量的关系，假设单价与批量的函数关系是一个多项式函数，可用多项回归来建立模型。7第7页/共93页6.1.4 模型的建立 根据模型的分析，可以建立多项式模型 令 ,则回归方程可写成 这是一个二元线性回归模型。且 ，其中：8第8页/共93页9第9页/共93页6.2 线性方程组迭代法概述迭代法概述迭代法概述10迭代法迭代法迭代矩阵迭代矩阵第10页/共93页迭代法：即用某种极限过程逐步逼近线性方程组精确解的方法。迭代法具有需要计算机存储较少、程序设计简单、原始系数矩阵在计算过程中始终不变等优点，但有收敛性或收敛速度的问题。迭代法是解大型稀疏矩阵方程组得重要方法。迭代法则能保持矩阵的稀疏性，具有计算简单，编制程序容易的优点，并在许多情况下收敛较快。故能有效地解一些高阶方程组。11第11页/共93页迭代法的基本思想 迭代法的基本思想是构造一串收敛到解的序列，即建立一种从已有近似解计算新的近似解的规则。由不同的计算规则得到不同的迭代法。12第12页/共93页 对线性方程组 ，其中 非奇异矩阵，。构造其形如 的同解方程组,其中M M为n n阶方阵，。13第13页/共93页 任取初始向量 ，代入迭代公式 产生向量序列 ，当k k充分大时，作为方程组AX=b AX=b 的近似解，这就是求解线性方程组的单步定常线性迭代法。M M 称为迭代矩阵。14第14页/共93页6.3 迭代法6.3.1Jacobi迭代法6.3.2Gauss-Seidel迭代法6.3.3迭代法的收敛性迭代法156.3.4超松弛代法第15页/共93页6.3.1 Jacobi迭代法对Ax=b Ax=b，设det(A)det(A)0 0，（6.3.1)6.3.1)将A A 改写成：(6.3.2)(6.3.2)16第16页/共93页 又将A A 分裂为：A=M-N A=M-N，则(6.3.1)(6.3.1)等价于 Mx=Nx+b(6.3.3)Mx=Nx+b(6.3.3)其中M M应选择为一个非奇异阵。并使Mz=f Mz=f 容易求解。对应(6.3.3)(6.3.3)可构造一个迭代过程：初始向量 ,(6.3.4)(6.3.4)特别地，若选取M=D M=D 则N=M-A=L+U N=M-A=L+U 从而(6.3.1)(6.3.1)化为：17第17页/共93页可得：JacobiJacobi迭代公式：（初始向量），其中：(6.3.5)(6.3.5)J J称为JacobiJacobi迭代的迭代矩阵。JacobiJacobi迭代的分量形式：引进记号：为第k k次近似，由6.3.5 6.3.5 有：18第18页/共93页19(6.3.6)(6.3.6)第19页/共93页 JacobiJacobi迭代公式简单，由公式(6.3.5),(6.3.6)(6.3.5),(6.3.6)可知，每迭代一次只需计算一次矩阵与向量乘法，计算机中只需要两组工作单元用来保存 及 且可用 来控制迭代终止。由迭代计算公式可知，迭代法一个重要特征是计算过程中原来矩阵 A A数据始终不变。20第20页/共93页例6.3.16.3.1 用JacobiJacobi迭代法求下面线形方程组，其精确 解是 ,21第21页/共93页解 先将Ax=b Ax=b 转化为等价方程组：迭代公式：选取初始向量 ,22第22页/共93页经1010次迭代解：误差为：23第23页/共93页 6.3.2 Gauss-Seidel迭代法 在(6.3.3)(6.3.3)中选取M=D-LM=D-L（下三角阵），则 N=M-A=U N=M-A=U，从而(6.3.1)(6.3.1)化为等价的(D-L)x(D-L)x=Ux+b (6.3.7)=Ux+b (6.3.7)可得Gauss-SeidelGauss-Seidel迭代公式：初始向量 ，其中 G G 称为Gauss-SeidelGauss-Seidel迭代矩阵24第24页/共93页G-SG-S迭代法的分量形式为：记有迭代公式(6.3.9)(6.3.9)，25第25页/共93页 Gauss-SeidelGauss-Seidel迭代法每迭代一次只需计算一次矩阵与向量的乘法，但G-SG-S迭代法比JacobiJacobi迭代法有一个明显的优点，那就是计算机上仅需一组工作单元用来保存分量 （或分量 ）当计算出 就冲掉旧的分量 。由G-SG-S迭代公式(6.3.9)(6.3.9)可看出在 的一步迭代中，计算分量 时利用了已计算出来的新分量 (j=1i1)(j=1i1)。因此，G-SG-S迭代法可以看作是JacobiJacobi迭代法的一个修正。26第26页/共93页例6.3.26.3.2 用G-SG-S方法解下面方程组，其精确解为：27第27页/共93页解 由(6.3.9)(6.3.9)可得本题G-SG-S迭代公式为：28第28页/共93页 经5 5次迭代得：从此例可见，G-S G-S迭代法比JacobiJacobi迭代法收敛快（初始向量相同，达到同样精度，所须迭代步数较少），但这个结论对Ax=bAx=b的矩阵 满足某些条件才是对的，甚至有这样的线性方程组，用JacobiJacobi方法是收敛的，而用G-SG-S迭代法却是发散的。29第29页/共93页 如线性方程组：此例用JacobiJacobi迭代法收敛而G-SG-S迭代法却发散。简要分析如下：30第30页/共93页31特征方程特征方程第31页/共93页6.3.3 迭代法的收敛性迭代法的迭代法的收敛性收敛性12条定理条定理8条定义32第32页/共93页 定义6.3.16.3.1 设有矩阵序列 及 如果 成立，则称 收敛于A A。记为 。例如，矩阵序列 当 时，（当 时）。矩阵序列收敛的概念可以用任何矩阵范数来描述。33第33页/共93页34下面介绍向量、矩阵的范数定义、常用的范数形式以及相关性质。定义 6.3.2 6.3.2 （向量范数）如果 的某个实值函数 （1 1）正定性：（2 2），c c 为实数（3 3）三角不等式：则称 为 上的一个向量范数（或向量的模）第34页/共93页 利用三角不等式容易证明：（4 4）在 中，三角不等式即为三角形两边长之和大于第三边长。如下图：35|x+y|y|第35页/共93页 定义 6.3.3 6.3.3设 上四种常用向量范数如下：(1)(1)向量的11范数：(2)(2)向量的 范数：(3)(3)向量的22范数：36第36页/共93页 (4)(4)向量的pp范数：易证明：它们均满足向量范数定义 ，且前三种向量范数均为第四种的向量范数的特例。其中 37第37页/共93页例6.3.3 6.3.3 设 计算 解 38第38页/共93页 定理 6.3.1 6.3.1 设非负实值函数 为 上任一向量范数，则N(x)N(x)是x x分量 的连续函数。定理 6.3.2 6.3.2（范数等价性）设 为 上任意两种范数，则存在常数 ，使得：定理 6.3.3 6.3.3 在 中，（时）其中 为向量的任一种范数39第39页/共93页定义 6.3.4 6.3.4 （矩阵的范数）若矩阵 的某个非负实值函数 满足条件：(1)(1)正定性：(2)(2)正齐次性：(3)(3)三角不等式：(4)(4)则称N(A)N(A)为 上的一个矩阵范数或模40第40页/共93页 如果将矩阵A A看成 向量空间中的元素，则由向量22范数定义有：这就是矩阵的FrobeniusFrobenius范数,明显它满足上面定义。在大多数与估算有关的问题中，矩阵和向量会同时参与讨论，所以希望引进一种矩阵的范数。它和向量范数相联系且和向量范数是相容的，即：41第41页/共93页 定义6.3.56.3.5 （矩阵的算子范数）设 给出一种向量范数 相应地定义一种矩阵的非负函数：（最大比值）易验证：满足前一矩阵范数的定义。因此 是 上的一个范数，称之为A A的算子范数。42第42页/共93页定理 6.3.4 6.3.4 设 是 上的一个向量范数，则 是 上的矩阵范数，且满足相容条件：43第43页/共93页 定理 6.3.5 6.3.5 设 则：（称为A A的行范数）(称为A A的列范数）（称为A A的22范数）其中 表示 的最大特征值。在计算上,（1 1），（2 2）比较容易，而（3 3）比较困难，但（3 3）在理论分析上十分有用。44第44页/共93页 定理 6.3.6 6.3.6 充要条件是 我们要考虑迭代法收敛性，即要研究B B在什么条件下使误差向量趋于零向量。即 定理 6.3.7 6.3.7 设 ，则 （零矩阵）的充要条件是：45第45页/共93页 定理 6.3.8 6.3.8 （迭代法基本定理）设有方程组 对于任意初始向量 及任意 f f，解此方程 组的迭代法（即 ）收敛的充要条件是 此定理应用时要计算谱半径，不太方便。将其修改成如下的形式：46第46页/共93页 定理6.3.96.3.9 （迭代法收敛的充分条件）设有方程组x=Bx+f x=Bx+f，且 为由迭代法产生的序列 ,为初始任意向量。若迭代矩阵B B 的某种范数满足 ，则：(1)(1)收敛于方程组 的唯一解 。(2)(2)（3 3）误差估计：47第47页/共93页 证明(1)(1)由于 ，即IB IB 可逆，有唯一解。设为 ，。引进误差向量 即得误差 的递推公式：于是 时 。48第48页/共93页由迭代公式有 又 则 即 利用此递推式即可得误差估计。49第49页/共93页例6.3.4 6.3.4 考察用JacobiJacobi方法求解线性方程组 的收敛性。50第50页/共93页解 则JacobiJacobi迭代矩阵为：51第51页/共93页 故解此方程组的JacobiJacobi方法收敛。52第52页/共93页 下面考察迭代法的收敛速度。，设B B 有n n 个线性无关的特征向量 ，相应的特征值为 ，由 （线性表示，基）可以看出当 愈小时，。愈快，即 愈快，故可用 来刻划迭代法的收敛快慢。53第53页/共93页现在来确定迭代次数k k，使 取对数得：定义6.3.66.3.6 称 为迭代法的收敛速度 由此可见 愈小，愈大，则所需的迭代次数就越少。54第54页/共93页 而n n 较大时，矩阵特征值计算比较困难，基本定理得条件比较难验证，所以最好建立与矩阵元素直接有关的条件来判断迭代法的收敛性。由于 ，所以可用 作为 上届的一种估计，这样的结果即为迭代法收敛的充分性条件。（如上面的定理）由这个充分性条件的结果(3)(3)可见 越小，收敛越快。另外，由其结果(2)(2)可知：55第55页/共93页 若 则 。因此可以用 作为控制迭代终止的条件。但要注意当 时，较大，尽管 很小，却有可能误差向量模 很大，迭代法收敛很慢而且此充分性条件的结果(3)(3)可以用来事先确定需要迭代多少次才能保证 。56第56页/共93页定理6.3.106.3.10 解方程组Ax=b Ax=b 的G-SG-S迭代法收敛的充分必要条 件是 ，G G为G-SG-S迭代矩阵。定义 6.3.7 6.3.7 （对角占优矩阵）设 若 ，即A A得每一行对角元素绝对值严格大于同行其他元素绝对值之和。则称A A为严格对角占优矩阵。若 且57第57页/共93页至少有一个不等式严格成立，则称 为弱对角占优矩阵。定义6.3.86.3.8 （可约与不可约矩阵）设 当 时，如果存在 n n阶排列矩阵p p 使 成立。其中 为 r r阶子矩阵，为阶子矩（），则称矩阵 A A为可约的。若不存在排列矩阵 p p使上式成立，则称 A A为不可约矩阵。58第58页/共93页 当A A 为可约矩阵，则 Ax=b Ax=b可经过若干行列重排化为两个低阶方程组求解。事实上，由 Ax=b Ax=b化为：，设 其中 为 r r 维向量。于是，求解Ax=b Ax=b 化为求解：另外A A 为可约矩阵的充要条件：存在指标集使 59第59页/共93页 定理6.3.116.3.11 （对角占优定理）：若 为严格对角占优阵或为不可约弱对角占优阵，则 A A是非奇异阵。证明(1)A(1)A为严格对角占优阵，采用反证法。若 det det A=0A=0，则Ax=0 Ax=0 有非零的解，设为 设 由奇次方程中的第k k 个方程得：则可得：60第60页/共93页即有：这与严格对角占优阵矛盾，故 61第61页/共93页(2)A2)A为不可约弱对角占优阵，采用反证法。设 使 Ax=0 Ax=0。并令m m 使 （由弱对角占优定义）成立。再定义下标集合：62第62页/共93页 在(6.3.10)(6.3.10)中取k=m k=m，将导致与(6.3.11)(6.3.11)得矛盾故 （空集合）。对任意 ，有 （由(6.3.10)(6.3.10)），由此可见，当 时有 。否 则，上式就与 A A为弱对角占优阵矛盾。但对任意 和 ，必有 ，因而 从而 A A为可约矩阵，这与A A 为不可约矩阵假设矛盾。63第63页/共93页定理 6.3.12 6.3.12 若 为严格对角占优矩阵或为不可约对角占优矩阵，则对任意的初始向量 方程组Ax=b Ax=b 的JacobiJacobi迭代法和G-SG-S迭代法均收敛。且G-SG-S迭代法比JacobiJacobi迭代法收敛速度快。证明 设A A为不可约对角占优矩阵，先证明G-SG-S迭代法收敛。由弱对角占优阵假设知 而G-SG-S迭代矩阵为 又 64第64页/共93页 65即即记记则：则：第65页/共93页 下面来证明当 时，这是因为A A 为不可约阵，则 C C也不可约，由 A A为弱对角矩阵，可得：且至少有一个不可约等式严格成立，这表明当 时，C C为不可约弱对角占优矩阵，于是：66第66页/共93页 由前一个定理可知当 时，这一结论表明det C=0 det C=0 的根 满足：即 故G-SG-S法收敛。在同一条件下，对于JacobiJacobi方法（）完全类似于上可证。当 A A为严格对角占优阵时，证明方法完全类似。67第67页/共93页6.3.4 超松弛代法 逐次超松弛迭代法（Successive Over Relaxation Successive Over Relaxation Method.Method.简称为SORSOR法）是Gauss-SeidelGauss-Seidel法的一种加速方法。这是解大型矩阵方程组的有效方法之一，具有计算公式简单，程序设计容易，占用计算机内存较少等优点，但需要选择好的加速因子，即最佳的加速因子。对 Ax=b (6.3.12)Ax=b (6.3.12)其中 为奇异矩阵，且设 分解：设已知第 k k次迭代向量 ，及第 k+1 k+1次迭代向量 的分量 现在来计算分量：68第68页/共93页先用G-SG-S迭代法求出辅助量 （预测）再取为与的某种平均值（即加权平均），即将(6.3.14)(6.3.14)代入(6.3.15)(6.3.15)即得解Ax=b Ax=b 的逐次超松弛迭代公式：69第69页/共93页其中称 为松弛因子，或写为：70第70页/共93页 显然 时，解(6.3.12)(6.3.12)的SORSOR法即为Gauss-Gauss-SeidelSeidel迭代法。SOR SOR法中每迭代一次，主要计算量是计算一次矩阵与向量乘法。由(6.3.16)(6.3.16)可见在计算机上用SORSOR法解方程组只需一组工作单位元，以便存放近似解。迭代算时，可用 来控制迭代终止。71第71页/共93页 当 时，称(6.3.16)(6.3.16)为低松弛法；当 时称(6.3.16)(6.3.16)为超松弛法。例 用SORSOR法解：其精确解为72第72页/共93页解 取 ，则SORSOR迭代公式为：取 ，第1111次迭代结果为：73第73页/共93页 对 取其它值，迭代次数如下表，由此例可见，松弛因子选择得好，会使SORSOR迭代的收敛大大加速。本例中 是最佳松弛因子。表6.3.1 6.3.1 结果表：74第74页/共93页75第75页/共93页下面考察SORSOR迭代公式的矩阵形式，由(6.3.16)(6.3.16)可改写为：由 ，则：即76第76页/共93页 由于对任意 均为奇异阵（设 而 L L为下三角阵，且对角线元素为0 0）则：因此，若 ，则SORSOR迭代公式为：其中 77（6.3.86.3.8）第77页/共93页 称 为SORSOR方法的迭代矩阵，应用关于迭代法的收敛性定理于(6.3.18)(6.3.18)可得：定理6.3.136.3.13 对 Ax=b Ax=b，且 则解方程组的充要条件是：。引进超松弛法的想法是希望能选择松弛因子使迭代过程(6.3.16)(6.3.16)收敛较快，也就是应选择因子 使 78第78页/共93页 下面考虑对于方程组 （6.3.126.3.12）（）超松弛因子在什么范围内取值才可能收敛。定理 6.3.14 6.3.14 （必要条件）对Ax=b ()Ax=b ()的SORSOR方法若收敛，则：。证明 由SORSOR法收敛及上定理，知 ，设 的特征值为 则：79第79页/共93页 80 而而因此因此 即即 。此结果表明要使此结果表明要使SORSOR法收敛法收敛，松，松弛因子弛因子 必须（必须（0 0，2 2）中。那么，反过来，若选取中。那么，反过来，若选取 在（在（0 0，2 2）中，）中，SORSOR法是法是否一定收敛呢？否一定收敛呢？即即第80页/共93页定理 （充分条件）若 A A为对称正定矩阵，且 则解（6.3.126.3.12）的SORSOR法必收敛。证明 若能证明在定理条件下，对 的任一特征值 有：，则定理得证。事实上，设y y为对应 的 的 特征向量。即：即 81第81页/共93页亦即：为找 的表达式，考虑内积：则有：82第82页/共93页 而：设 ，由于 则 ，从而：83第83页/共93页因此：从而：当 由 （6.3.206.3.20），（6.3.216.3.21）有：从而由(6.3.21)(6.3.21)可知 ，故SORSOR收敛。84第84页/共93页 6.4 关于回归模型的求解在6.16.1的问题中，令，85第85页/共93页则可以得到一个线性方程组：令86第86页/共93页用Gauss-SeidelGauss-Seidel迭代法87第87页/共93页 用MATLABMATLAB编写Gauss-SeidelGauss-Seidel迭代法的程序，经过50005000次迭代得到结果如下：所以二次回归模型为：88第88页/共93页习题 61 1、取初始向量X(0)=(0,0,0)T,X(0)=(0,0,0)T,用雅可比迭代法求解线性方程组89第89页/共93页习题 62 2、用Gauss-SeidelGauss-Seidel迭代法解线性方程组：90第90页/共93页习题 63 3、取 ，用SORSOR迭代法解线性方程组：91第91页/共93页中南大学数学科学与计算技术学院中南大学数学科学与计算技术学院第92页/共93页感谢您的观看！第93页/共93页