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

《数值线性代数课程设计》专业：信息与计算科学班级： 13405011学号： **********姓名：***实验日期： 2016.05.09报告日期： 2015.05.13实验地点：数理学院五楼机房超定方程组的求解邢耀光（班级：13405011 学号1340501123）摘要：在实验数据处理和曲线拟合问题中，求解超定方程组非常普遍。

比较常用的方法是最小二乘法。

形象的说，就是在无法完全满足给定条件的情况下，求一个最接近的解。

最小二乘法（又称最小平方法）是一种数学优化技术。

它通过最小化误差的平方和寻找数据的最佳函数匹配。

利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。

关键字：最小二乘问题，残量，超定方程组，正则化方程组，Cholesky 分解定理。

正文：最小二乘法的背景：最小二乘法（又称最小平方法）是一种数学优化技术。

它通过最小化误差的平方和寻找数据的最佳函数匹配。

利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。

最小二乘法还可用于曲线拟合。

其他一些优化问题也可通过最小化能量或最大化熵用最小二乘法来表达。

最小二乘法经常运用在交通运输学中。

交通发生预测的目的是建立分区产生的交通量与分区土地利用、社会经济特征等变量之间的定量关系，推算规划年各分区所产生的交通量。

因为一次出行有两个端点，所以我们要分别分析一个区生成的交通和吸引的交通。

最小二乘问题：最小二乘问题多产生于数据拟合问题。

例如，假定给出m 个点1,...,m t t 和这m 个点上的实验或观测数据1,...,my y ，并假定给出在i t 上取值的n 个已知函数1(),...,()n t t ψψ。

考虑i ψ 的线性组合1122(;)()()...()n n f x t x t x t x t ψψψ=+++ ， （1）我们希望在1,...,m t t 点上(;)f x t 能最佳的逼近1,...,m y y 这些数据。

为此，若定义残量1()()ni i j j i j r x y x t ψ==-∑ ， 1,...,i m = ， （2）则问题成为：估计参数1,...,n x x ，使残量1,...,m r r 尽可能地小。

（2）式可用矩阵-向量形式表示为()r x b Ax =- ， （3） 其中1111()(),()()n m n m t t A t t ψψψψ⎛⎫ ⎪= ⎪ ⎪⎝⎭ 1,m y b y ⎛⎫⎪= ⎪ ⎪⎝⎭1(,...,),T n x x x = 1()((),...,()).T m r x r x r x = 当m n =时，我们可以要求()0r x =，则估计x 的问题就可以用第一章中讨论的方法解决。

当m n >时，一般不可能使所有残量为零，但我们可要求残向量()r x 在某种范数意义下最小。

最小二乘问题就是求x 使残向量()r x 在2范数意义下最小。

定义1：给定矩阵m n A R ⨯∈及向量m b R ∈，确定n x R ∈，使得2222()min ()min .nny R y R b Axr x r y Ay b ∈∈-===- （4）这就是所谓的最小二乘问题，简称为LS 问题，其中的()r x 常常被称为残向量。

在所讨论的最小二乘问题中，若r 线性依赖于x ，则称其为线性最小二乘问题：若r 非线性依赖于x ，则称其为非线性最小二乘问题。

最小二乘问题的解x 又可称做线性方程组,Ax b = m n A R ⨯∈ （5） 的最小二乘解，即x 在残向量()r x b Ax =-的2范数最小的意义下满足方程组（5）。

当m n >时称（5）式为超定方程组。

定理1:（Cholesky 分解定理） 若n n A R ⨯∈对称正定，则存在一个对角元均为正数的下三角阵n n L R ⨯∈，使得.TA LL = （6） （6）式称为Cholesky 分解，其中的L 称作A 的Cholesky 因子。

因此，若线性方程组Ax b =的系数矩阵是对称正定的，则我们自然可按如下的步骤求其解： （1）计算A 的Cholesky 分解：T A LL = ； （2）求解Ly b =得y ；（3）求解TL x y =得x ； 简单而实用的方法是直接比较T A LL =两边的对应元素来计算L 。

设11212212.n n nn l l l L l l l ⎛⎫⎪⎪= ⎪⎪⎝⎭比较TA LL =两边对应的元素，得关系式1jij ip jpp a ll ==∑ ，1j i n ≤≤≤ （7）首先，由21111a l =，得11l = 再由1111i i a l l =，得1111i i l a l =，1,...,.i n =这样便得到了矩阵L 的第一列元素。

假定已经算出L 的前1k -列元素，由 21,kkk kpp a l==∑得11221.k kk kk kp p l a l -=⎛⎫=- ⎪⎝⎭∑ （8）再由11k ik ip kp ik kk p a l l l l -==+∑ ， 1,...,,i k n =+得11k ik ik ip kp kk p l a l l l -=⎛⎫=- ⎪⎝⎭∑ ，1,...,.i k n =+ （9）这样便求出了L 的第k 列元素。

这种方法称为平方根法。

记最小二乘解的解集为LS χ，即}{:n LS x R x LS χ=∈是问题(3)的解 ，定理2：LS x χ∈ 当且仅当.T T A Ax A b = （10）方程组（10）常常被称为最小二乘问题的正则化方程组或法方程组，它是一个含有n 个变量和n 个方程的线性方程组。

在A 的列向量线性无关的条件下，T A A 对称正定，故可用平方根法求解方程组（6），这样，我们就得到了求解最小二乘问题最古老的算法———正则化方法，其基本步骤如下： （1）计算,TC A A = ;Td A b =（2）用平方根法计算C 的Cholesky 分解：;TC LL = （3）求解三角方程组Ly d =和.TL x y =实验 ：一：超定方程组的求解原理：设A 是m n ⨯阶矩阵()m n >，则线性方程组Ax b =为超定方程组，这里,mmx R b R ∈∈。

如果A 的秩为n ，则称A 为列满秩矩阵。

超定方程组的解满足法方程T TA Ax A b =，该解使得22min b Ax-=，称之为最小二乘解。

题目： 2222211 1.1 1.1211.2 1.2311.3 1.341 1.4 1.451 1.51.5x ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦用正则化方法求解，要求:（1）T B LL = 不得使用MathCAD Cholesky 指令； （2）TB LL =使用MathCAD Cholesky 指令。

解：（1）222221 1.1 1.111.2 1.211.3 1.31 1.4 1.41 1.51.5A ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦, 则 5 6.58.556.58.5511.3758.5511.37515.298TB A A ⎡⎤⎢⎥==⎢⎥⎢⎥⎣⎦1520.528.25T g A b ⎡⎤⎢⎥==⎢⎥⎢⎥⎣⎦,11 2.236L ==, 2121112.907B L L == ,3131113.824B L L == ,220.316L == ,32312132220.822B L L L L -== ,330.037L == ,即 2.236002.9070.31603.8240.8220.037L ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦ , 2.236 2.907 3.82400.3160.822000.037TL ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦,则1136.7083.1629.27310y L g --⎡⎤⎢⎥==⎢⎥⎢⎥⨯⎣⎦ , ()11110102.47810Tx L y ---⎡⎤⎢⎥==⎢⎥⎢⎥⨯⎣⎦,x 即为所求的最小二乘解。

（2）222221 1.1 1.111.2 1.211.3 1.31 1.4 1.41 1.51.5A ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦， 2.23600() 2.9070.31603.8240.8220.037cholesky B ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦ ,则 2.236002.9070.31603.8240.8220.037L ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦， 2.236 2.907 3.82400.3160.822000.037TL ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦,则1136.7083.1629.27310y L g --⎡⎤⎢⎥==⎢⎥⎢⎥⨯⎣⎦ , ()11110102.47810Tx L y ---⎡⎤⎢⎥==⎢⎥⎢⎥⨯⎣⎦,x 即为所求的最小二乘解。

二：已知如下数据：12.y a x a =+利用最小二乘法拟合曲线解：令0.00.20.40.60.8 1.0 1.20.9 1.9 2.8 3.3 4.0 5.7 6.5B ⎡⎤=⎢⎥⎣⎦ ,0.00.20.40.60.81.01.2x ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦ ，0.91.92.83.34.05.76.5y ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦则10.010.210.410.610.81 1.01 1.2A ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦,()10.8434.571TT X A A A y -⎡⎤==⎢⎥⎣⎦ ，即()0.843 4.571p x x =+ ， 故最小二乘法拟合曲线为 4.5710.843.y x =+程序附录： 一；,, ,,, , ,,A 111111.11.21.31.41.51.121.221.321.421.52⎛ ⎝⎫⎪⎪⎪⎪⎪⎪⎪⎪⎭:=b 12345⎛ ⎝⎫⎪⎪⎪⎪⎪⎭:=Ax b:=B A T A :=g A T b :=B 56.58.55 6.58.5511.3758.5511.37515.298⎛ ⎝⎫⎪⎪⎭=g 1520.528.25⎛ ⎝⎫⎪⎪⎭=f B ()n rows B ()←L identity n ()←Lk k,Bk k,←k1if Lk k,Bk k,1k 1-p L k p,()2∑=-←otherwiseLi k,Bi k,Lk k,←k 1if break ()knif Li k,Bi k,1k 1-p L i p ,L k p ,()∑=-Lk k,←otherwisei k 1+n..∈for k 1n ..∈for L:=,,,,, ,, ，二；, , ,, , , , ,, , ,,,, ,, .()0.843 4.571.p s x =+心得体会：通过本次的课程设计，让我学会了很多，学会了简单的MathCAD 软件的用法。