数值线性代数课程设计—超定方程组的求解

1、数值线性代数课程设计专业：信息与计算科学班级：13405011学号：1340501123姓名：邢耀光实验日期：2016.05.09报告日期：2015.05.13实验地点：数理学院五楼机房超定方程组的求解邢耀光(班级：13405011 学号 134050112摘要：在实验数据处理和曲线拟合问题中，求解超定方程组非常普遍。 比较常用的方法是最小二乘法。形象的说，就是在无法完全满足给定条件的情况下，求一个最接近的解。最小二乘法(又称最小平方法) 是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地 求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和

2、为最小。关键字： 最小二乘问题，残量，超定方程组，正则化方程组，Cholesky分解定理。正文：最小二乘法的背景：最小二乘法(又称最小平方法)是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函 数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方 和为最小。最小二乘法还可用于曲线拟合。其他一些优化问题也可通过最小化能量或最大化嫡用最小二乘 法来表达。最小二乘法经常运用在交通运输学中。交通发生预测的目的是建立分区产生的交通量与分区土地利用、社会经济特征等变量之间的定量关系，推 算规划年各分区所产生的交通量。因为一次出行有两个端点，所以我们要分别分

3、析一个区生成的交通和吸 引的交通。最小二乘问题：最小二乘问题多产生于数据拟合问题。例如，假定给出m个点t1,.,tm和这m个点上的实验或观测数据V，ym 并假定给出在t上取值的n个已知函数 H,.Wn(t)。考虑 的线性组合f (x;t) =Xi巴 +X2中 2(t) + + n(t) ,(1)我们希望在".,。点上f(x;t)能最佳的逼近 V,.,ym这些数据。为此，若定义 残量nr,(x) =y, -Z xj j(ti) , i=1,.,m ,(2)j 1则问题成为：估计参数 为,.，4,使残量1,.,如尽可能地小。(2)式可用矩阵-向量形式表示为其中(3)r(x) = b -

4、Ax ,俚1也)%(L)'A 二 |：M(tm)%(tm)/x = (x，. XT )r(x) =(r1(x),.,rm(x)T.当m = n时，我们可以要求r(x)=0,则估计x的问题就可以用第一章中讨论的方法解决。当 mn时,般不可能使所有残量为零，但我们可要求残向量 r(x)在某种范数意义下最小。 最小二乘问题就是求 x使残向量r(x)在2范数意义下最小。定义1:给定矩阵A乏Rm"n及向量bw Rm ,确定xw Rn ,使得b-Ax2 . r(x)2 .rmin=愧 Ay - b 2.y ny n这就是所谓的 最小二乘问题，简称为LS问题，其中的r(x)常常被称为残向量

5、。(4)在所讨论的最小二乘问题中，若r线性依赖于x ,则称其为线性最小二乘问题：若r非线性依赖于x,则称其为非线性最小二乘问题。最小二乘问题的解x又可称做线性方程组Ax=b, AwRm 刈(5)的最小二乘解，即 x在残向量r(x) =bAx的2范数最小的意义下满足方程组(5)。当mn时称(5) 式为超定方程组。定理1： ( Cholesky分解定理) 若AwRn对称正定，则存在一个对角元均为正数的下三角阵L w Rn* ,使得A =LLT.(6)式称为Cholesky分解，其中的L称作A的Cholesky因子。(6)因此，若线性方程组 Ax =b的系数矩阵是对称正定的，则我们自然可按如下的步骤

6、求其解:(1)计算 A 的 Cholesky 分解：A = LLT ；(2)求解Ly =b得y ；(3)求解 LT x = y 得 x ；简单而实用的方法是直接比较A = LLT两边的对应元素来计算,In、l21 l22L = :l /l c lI n11 n25n ,比较A = LLT两边对应的元素，得关系式aij . l ipl jpp=1首先，由a11 =l112,得l11a11.再由ai1 T11L1 ，得li1 = ai1, l11 , i - 1>.1 n.这样便得到了矩阵L的第一列元素。假定已经算出L的前k -1列元素，由kakk = '、lkp ，P 1得1k&q

7、uot;2lkk akk 一乙 1kp .1Rd1再由k 1dkliplkp+liklkk , i=k+1,.,n,p 1得<k)/l ik aik 乙 l ipl kp /lkk , i k 1,., n.1pW)1(8)(9)这样便求出了 L的第k列元素。这种方法称为 平方根法。记最小二乘解的解集为 ？LS,即工LS =xw Rn ： x是LS问题(3)的解,定理2: xwls当且仅当A Ax = A b.(10)方程组(10)常常被称为最小二乘问题的 正则化方程组或法方程组，它是一个含有n个变量和n个方 程的线性方程组。在 A的列向量线性无关的条件下，atA对称正定，故可用平方根法

8、求解方程组(6),这样，我们就得到了求解最小二乘问题最古老的算法正则化方法，其基本步骤如下：(1)计算 C =AtA, d =ATb;(2)用平方根法计算 C的Cholesky分解：C = LLT;(3)求解三角方程组 Ly =d和LTx = y.实验：一：超定方程组的求解原理：设A是mxn阶矩阵(mn),则线性方程组 Ax=b为超定方程组，这里 xw Rm,bw Rm。如 果A的秩为n ,则称A为列满秩矩阵。超定方程组的解满足法方程AT Ax = ATb ,该解使得2|b -Ax|2 =min ,称之为最小二乘解。题目:一1 1.1卜1.21 1.31 1.4：1 1.51.111.221.

9、32 x1.421.52一112345J用正则化方法求解，要求：(1) B =LLt 不得使用 MathCAD 指令 Cholesky ;(2)B = LLt 使用 MathCAD 指令 Cholesky。11解：（1） A = 111J1.1 1.121.2 1.221.3 1.321.4 1.421.5 1.52一 56.58.55 T则 B=AtA= 6.58.5511.375：8.5511.37515.2981- 15 g =ATb = | 20.5 , L11 =7B1 =2.236, L21 = B21 = 2.907 ,28.25 一L11Lii= 3.824L22=辰-U f

10、=0.316B32 - L31 L21L32 二1-22= 0.822 ,-33 B33 - -3122-32= 0.037 ,2.236即 L = 2.9073.82400.3160.8220.0376.7083.162了.273父10,3 _Lt2.23602.9070.3163.8240.8220.037-10x= LT，y =102.478x101 -x即为所求的最小二乘解。1.11.21.31.41.51.121.221.321.421.522.2362.907?824cholesky(B)0.3160.8220 00.0372.2362.9073.824Lt0.3160.822 0

11、.0372.23602.9070.3163.8240.8220.037,6.708-103.16210_ _ _13 9273 M10_11 工.478父10 _x即为所求的最小二乘解。二：已知如下数据:X0.00.20.40.60.81.01.2V0.91.92.83.34.05.76.5利用最小二乘法拟合曲线y = a1x + a2.0.0 0.2解：令B =0.9 1.9-0.01-0.910.21.90.42.80.4 0.6 0.8 1.0 1.21,x =0.6,y =3.32.8 3.3 4.0 5.7 6.5一0.84.01.05.71.2i i6.5_1 0.01 0.21

12、0.4则 A= 1 0.6 ,1 0.81 1.01 1.2故最小二乘法拟合曲线为T 1 T=ATAATy =0.843U.5711即 p(x) =0.843 4.571xy =4.571x 0.843.程序附录:11A ：= 121.1 1.11.2 1.221.3 1.32 | b21.4 1.4iT、2：=346.58.55121.5 1.5 )f(B)一n 一 row s(B)L - iden tityn )for k三1 ,nif k = 1kfor i 三k . 1 ,nLkk 1. iLkP =1(b reak)Li k -Ax 尸b B :=ATAg =ATb, B =6.58

13、.558.55 1 1.3 751 1.3 75 1 5.2 98,15，20.5128.25Bi k if k = 1Lk kif k - nBi k , 一 (LiLk k2,oth erw isePPLk Poth erw ise-2.23600力.23600ltf (B)=2.9070.3160L =2.9070.3160口.8240.8220.037 jL ：=f(B) ,8240.8220.037 j,-2.23 6 2.90 7 3.82 400.31 6 0.82 21 000.03 7；_10106.708 3.162y ;=L_1g , x二 LT -1y-139.273 10-11 2.478 x 10,x 二 BT0.0 0.2 0.4 0.6 0.8 1.0 1.2B =0.9 1.9 2.8 3.3 4.0 5.7 6.5<0、0.20.4x = 0.60.81I2,t.9 11.92.83.345.7©5尸y, A =111111x1 =0 , x2 =0.2 , n =rows(x) , n =7 , i =1 ,n , Ai 1 ：=100.20.40.60.811.2>-1,x 尸:!ataT0.843A y ,X, p(s)4.571p(s) =0.843 4.571x