数值分析(最小二乘法)模板课件.ppt

1、1/4614:50最小二乘解的存在唯一性最小二乘解的存在唯一性最小二乘解的数值方法最小二乘解的数值方法数值分析 162/4614:50 x x1 x2 xm f(x)y1 y2 ym离散数据的直线拟合离散数据的直线拟合求拟合函数求拟合函数:mmyyyccxxx212121111 Ac=y超定方程组超定方程组xccx21)(1121yxcc 2221yxcc mmyxcc 213/4614:50 x x1 x2 xm f(x)y1 y2 ym离散数据的多项式拟合离散数据的多项式拟合求拟合函数求拟合函数:1011211nnmmnmyaxxyxxay 01()nnxaa xa x 超定方程组超定方程

2、组4/4614:50 x x1 x2 xm f(x)y1 y2 ym离散数据的线性拟合离散数据的线性拟合求拟合函数求拟合函数:1001111201()()()()()()nmmnmnmyaxxxyxxxay 0011()()()()nnxaxaxax超定方程组超定方程组5/4614:50,m nAxbARmn 超超定定方方程程其其中中 回顾回顾:22 argmin|xxAnextbesbtx 最最小小二二乘乘解解（）22argmin()TTxAxbA AxA bnormal eqaution初初等等变变分分原原理理 1,()TTTA AxA AA b 进进一一步步地地如如果果可可逆逆 则则 A

3、xb 6/4614:50最小二乘拟合问题研究包括最小二乘拟合问题研究包括:模型的选取模型的选取存在唯一性存在唯一性最小二乘解的计算最小二乘解的计算7/4614:50广义矩阵广义矩阵(Ax=b统一的理论解释统一的理论解释)Axb=Ax b8/4614:50相容方程的解相容方程的解定义定义:一个方程组称为相容方程一个方程组称为相容方程(consistent equation),若若至少存在一个解能够严格满足该方程组。至少存在一个解能够严格满足该方程组。定理定理:线性方程线性方程Ax=b是相容的当且仅当增广矩阵的秩是相容的当且仅当增广矩阵的秩等于矩阵等于矩阵A的秩的秩,即即rank(A,b)=ran

4、k(A)。定理定理:相容方程相容方程Ax=b对对y不等于零有解不等于零有解x=Gb当且仅当当且仅当AGA=A。(G称为是称为是A的广义逆的广义逆generalized inverse)-1-1=,()min,),m nAxA bmnArank Am nxA bxGb 如如果果且且 非非奇奇异异 则则方方程程的的解解为为。一一个个自自然然的的问问题题是是在在和和 为为秩秩亏亏缺缺的的情情况况下下是是否否存存在在一一个个与与相相类类似似的的解解 比比如如是是相相容容方方程程的的解解?9/4614:50相容方程解的唯一性相容方程解的唯一性是否存在某种意义下的唯一性是否存在某种意义下的唯一性?22=m

5、in,Ax bGGbxGbG 如如果果存存在在 满满足足则则称称为为相相容容方方程程的的最最小小范范数数解解广广义义逆逆矩矩阵阵 为为最最小小范范数数广广义义逆逆矩矩阵阵。最小范数解最小范数解(minimum norm solution):定理定理:Gb是相容矩阵的最小范数解当且仅当是相容矩阵的最小范数解当且仅当 AGA=A,(GA)H=GA。参考参考:张贤达张贤达,矩阵分析与应用矩阵分析与应用,清华大学清华大学10/4614:50不相容方程解的存在性不相容方程解的存在性不相容方程的最小二乘解总是存在的。不相容方程的最小二乘解总是存在的。证明证明:即证明正规方程是相容方程。即证明正规方程是相容

6、方程。rank(ATA,b)=rank(ATA)rank(),rank()rank()rank(),rank(,)rank(),rank(,)minrank(),rank(,)rank(,)k,rank(,)rank()TTTTTTTTTTTTTTTTAkA AAAkA A A bA AkA A A bAA bA A A bAA A A bA bA AkA bA A 设设则则由由于于故故综综上上所所述述Axb 22 argmin|xxAxb11/4614:50定理定理 如果矩阵如果矩阵A 列满秩列满秩,则则ATA可逆可逆。121122:,0,0,nnnTTTTcAccccc A AcA AA

7、A 如如果果矩矩阵阵列列满满秩秩则则矩矩阵阵列列向向量量线线性性无无关关 则则对对于于任任意意的的非非零零向向量量进进一一步步有有对对任任意意非非零零向向量量因因为为矩矩阵阵正正定定证证明明可可逆逆。定理定理 矩阵矩阵A列满秩时列满秩时,最小二乘解唯一最小二乘解唯一x=(ATA)-1ATb。12/4614:50不相容方程解的唯一性不相容方程解的唯一性是否存在某种意义下的唯一性是否存在某种意义下的唯一性?2222G:,GbxxxAxbAzbzGbG若若存存在在 满满足足其其中中则则称称是是最最小小范范数数最最小小二二乘乘解解称称为为最最小小范范数数最最小小二二乘乘广广义义矩矩阵阵。最小范数最小二

8、乘解最小范数最小二乘解(minimum norm least squares solution)定理定理:Gb是不相容矩阵的最小范数最小二乘解当且仅当是不相容矩阵的最小范数最小二乘解当且仅当 AGA=A,(AG)H=AG,GAG=G,(GA)H=GA。注释注释:最小范数最小二乘广义矩阵即最小范数最小二乘广义矩阵即Moore-Penrose矩阵。矩阵。13/4614:50总结总结相容方程相容方程 矩阵可逆则解唯一矩阵可逆则解唯一,如果矩阵秩亏损的情形如果矩阵秩亏损的情形,则所有解则所有解中有唯一的最小范数解。中有唯一的最小范数解。不相容方程不相容方程 首先最小二乘解一定存在首先最小二乘解一定存在

9、,如果矩阵列满秩则最小二如果矩阵列满秩则最小二乘解唯一乘解唯一,如果矩阵秩亏损的情形如果矩阵秩亏损的情形,所有最小二乘解有所有最小二乘解有唯一的最小范数最小二乘解。唯一的最小范数最小二乘解。Axb Axb 22argminxAxb 14/4614:50对于任意矩阵对于任意矩阵,Moore-Penrose逆矩阵存在且唯一。逆矩阵存在且唯一。Matlab:pinv(Pseudoinverse)111,(),(),TTTTXXXXXX XXXXXXX 如如果果 是是方方阵阵且且非非奇奇异异 则则如如果果 是是列列满满秩秩如如果果 是是行行满满秩秩比较比较back slash和和pinv的区别。的区别

10、。1231645617,789181011121913141520XyXy,pinv(X)*y,norm(Xy),norm(pinv(X)*y)15/4614:50参考文献参考文献:Sparse and Redundant Representations:From Theory to Applications in Signal and Image Processing1范范数数意意义义下下的的残残差差最最小小16/4614:50最小二乘拟合问题研究包括最小二乘拟合问题研究包括:模型的选取模型的选取存在唯一性存在唯一性最小二乘解的计算最小二乘解的计算17/4614:50为什么不直接求解正规方程

11、为什么不直接求解正规方程?()TTnormal equAationAxA b 正正规规方方程程 ,m nTn nTTA AxA bARmnA AR 其其中中注注意意2()()Tcond A Acond A 1011211nnmmnmyaxxyxxay 2211110,110TXX X 18/4614:50初等行变换不改变方程组的解初等行变换不改变方程组的解1.交换矩阵第交换矩阵第i行与第行与第j行行2.非零数非零数k乘以矩阵第乘以矩阵第i行的每个元素行的每个元素3.矩阵第矩阵第i行的每个元素的行的每个元素的k倍加到第倍加到第 j行的对应元素行的对应元素12113112111xx 1233911

12、2111xx 1.7500 0.7500 1.9500 0.950019/4614:50A(n 1)=Fn-1Fn-2F1 A其中其中Fk 为为 Frobenius矩阵。矩阵。A=F1-1F2-1 Fn-1-1 A(n 1)直接方法直接方法:高斯消元法高斯消元法 L U矩阵矩阵LU分解是高斯消元法的矩阵编码。分解是高斯消元法的矩阵编码。1111,121nnnmmmL )1()1(2)1(2211211nnnnnaaaaaaU20/4614:50,m nAxbARmn 超超定定方方程程其其中中 回顾回顾:22 argmin|squaresxxleAxbast最最小小二二乘乘解解（）22argmi

13、n()TTxAxbA AxA bnormal eqaution初初等等变变分分原原理理 Axb 不不相相容容21/4614:50回顾回顾:正交矩阵乘向量正交矩阵乘向量,则向量则向量2范数不变。范数不变。QTQ=I,y=Qx()()TTTTTy yQxQxx Q Qxx x222222|yQxx2222 argmin|argmin|xxRxQAxQbxQb正交矩阵正交矩阵QTQ=QQT=I半正交矩阵半正交矩阵QTQ=I(列正交列正交)或或QQT=I(行正交行正交)22/4614:50 Gram-Schmidt正交化正交化12111323112211111112122131311(,)(,)(,)

14、(,)(,)(,)(,-)(,)(,)(,)kkkkku vu uu vuvu uuuu vukuukuukvuuuuuuuuvvvvu 23/4614:50 Gram-Schmidt正交化正交化12111323112211111112112123113(,)(,)(,)(,)(-,)(,)(,)(,)(,)(,)+kkkkku vu uu vuvu uuuku vuvu uuukkuuuuuvvvuuuvu 24/4614:50Gram-Schmidt正交化的矩阵编码正交化的矩阵编码12111323112211111112112123113(,)(,)(,)(,)(-,)(,)(,)(,)(

15、,)(,)+kkkkku vu uu vuvu uuuku vuvu uuukkuuuuuvvvuuuvu 12131232121121111,kkkkkkrrrrv vvru uur u1,u2,un是正交基向量是正交基向量R单位上三角矩阵单位上三角矩阵25/4614:50 Gram-Schmidt正交化正交化121113231122111111121311112222212232233131(,)(,)(,)(,)(,)(,)(,)(,)(,)(,)-,/|,/|,/|,/|kkkkkq vq qq vqvq qqqq vqvq qqkkqkkkkquvvvuuquuquuqqququuq

16、uuvq 121113231122111111122232121312112213(,)(,)(,)(,)(,)(,)(,)(,)(,)-(,)|kkkkkq vq qq vqvq qqqq vqkkvkq qqqkquuuqqvqvqqqvquqv 26/4614:50 Gram-Schmidt正交化正交化121113231122111111112123112212122332(,)(,)(,)(,)(,)(,)(,)(,)(-,),)|=+|+|+kkkkkq vq qq vqvq qqqq vqvkq qqqkkkqqqqqqqquuuquvvvv 1212221323112123123

17、1311212-|(,)|(,)(,)|(,)(,)|+|=+kkkkkkkvvvqqqqquq vuq vq vuqqqvqqvvqu 27/4614:50Gram-Schmidt正交化的矩阵编码正交化的矩阵编码12121312223232121122,|kkkkkkkurrrurruq qqruv vv q1,q2,qn是标准正交基向量是标准正交基向量,Q正交矩阵正交矩阵,R上三角矩阵上三角矩阵Matlab命令命令:qr1122132331121211123311-=+|(,)|(,)(,)|(,)(,)|+|kkkkkkkuq vuq vq vqqqqqqquq vqvvvuvqqv 2

18、8/4614:50矩阵的正交三角分解矩阵的正交三角分解:A=QR,(,0),()m nn nm mm nAQRRRq rqr AAQRRRRq rqrRAQQ精精简简版版本本完完整整版版本本注释注释:经典的经典的Gram-Schmidt过程数值稳定过程数值稳定性不令人满意性不令人满意,一般将一系列一般将一系列Householder变换变换(正交变换正交变换)作用于最小二乘问题。作用于最小二乘问题。QR分解有两种版本分解有两种版本:完全完全QR分解和精简分解和精简QR分解。分解。29/4614:50A=2 4;3-5;1 2;Q,R=qr(A);%完整型完整型y=11 3 6;x=R(Q*y);

19、norm(A*x-y)21229.88873.74171.3363argmin 8.247106.5738 0.470042xx 243512A1136y 解法解法1(完整完整QR分解分解)2222argminargmin,TTTQ QRxQ yRxQ yQ其其中中 是是正正交交矩矩阵阵。30/4614:50243512A1136y 解法解法2(精简精简QR分解分解)Q,R=qr(A,0);%精简型精简型x=R(Q*y);norm(A*x-y)21223.74171.33639.8887argmin06.5738 8.2471xx 2211122122211222argminargminarg

20、min,TTTTTTTTQQRxQ yQ RxyRxQ yQ QQQQ y 其其中中 是是正正交交矩矩阵阵。31/4614:5021229.88873.74171.3363argmin 8.247106.5738 0.470042xx 完整和精简完整和精简QR分解的比较分解的比较21223.74171.33639.8887argmin06.5738 8.2471xx 2222argminargmin,QTTTQ QRxQ yRxQ y其其中中 是是正正交交矩矩阵阵。2211122122211222argminargminargmin,TTTTTTTTQQRxQ yQ RxyRxQ yQ QQQ

21、Q y 其其中中 是是正正交交矩矩阵阵。32/4614:50Ab (算法算法QR分解分解,具体实现具体实现Household变换变换)Tools Basic FittingIf A is an M-by-N matrix with M N or NM then X=AB is the solution in the least squares sense to the under-or overdetermined system of equations A*X=B.Matlab超定方程组求解超定方程组求解Matlab拟合拟合GUI33/4614:50Matlab的多项式拟合命令的多项式拟合命

22、令polyfit和和polyval调用格式调用格式P=polyfit(X,Y,N)调用格式调用格式Y=polyval(P,X)load census,plot(x,y,o)P=polyfit(x,y,2);polyval(P,2010)截至截至2010年年4月月1日日,美国居住人口总数为美国居住人口总数为308,745,53834/4614:50Matlab的非线性拟合命令的非线性拟合命令非线性拟合非线性拟合 lsqnonlin和和 lsqcurvefitx=0:.1:10;y=0.12*exp(-0.213*x)+0.54*exp(-0.17*x).*sin(1.23*x);f=inline

23、(a(1)*exp(-a(2)*x)+a(3)*exp(-a(4)*x).*sin(a(5)*x),a,x);a,res=lsqcurvefit(f,1,1,1,1,1,x,y)241135102argmin(sin()kkxxkkkaaaeexyaaa 35/4614:50严格满足插值条件严格满足插值条件 vs 追求最小残差平方和追求最小残差平方和适定方程组求解适定方程组求解 vs 最小二乘问题求解最小二乘问题求解寻找数据的规律寻找数据的规律(函数函数)或者说是压缩数据或者说是压缩数据思考思考:插值与拟合的异同插值与拟合的异同36/4614:50作业作业题目题目:数据的挖掘数据的挖掘(拟合或

24、插值拟合或插值)数据数据+挖掘挖掘知其然知其然,知其所以然知其所以然,用其然用其然,利其然利其然37/4614:50测试数据集测试数据集 Statistical Reference Dataset(http:/www.itl.nist.gov/div898/strd/)38/4614:50例例1 饮料的定价策略饮料的定价策略 一家公司在一家公司在22个近似相等大小的城市尝试销售一种新型个近似相等大小的城市尝试销售一种新型的运动型饮料的运动型饮料,售价售价(美元美元)以及在城市中每周的销量如下表以及在城市中每周的销量如下表:如果每件产品的制造成本是如果每件产品的制造成本是0.23美元美元,公司如

25、何设置公司如何设置全国统一售价全国统一售价 利润最大化利润最大化?城市城市售价售价销量销量/周周12345678910110.590.800.950.450.790.990.900.650.790.690.7939802200185061002100170020004200244033002300城市城市售价售价销量销量/周周12131415161718191021220.491.090.950.790.650.450.600.890.790.990.85600119019602760433069604160199028601920216039/4614:50例例2 非诚勿扰攻略非诚勿扰攻略例

26、例3 微信微信3点定位的真与假点定位的真与假 2012年年11月月4日日,一条微博称微信可以通过一条微博称微信可以通过三点定位法确定使用者的位置三点定位法确定使用者的位置,即记住自即记住自己的位置和与某人之间的距离己的位置和与某人之间的距离,变换两次变换两次位置重新记录距离位置重新记录距离,以这三个点为圆心、以这三个点为圆心、距离为半径画圆距离为半径画圆,交点就是要找的人的位交点就是要找的人的位置置,圆圈越多圆圈越多,位置越精确位置越精确。提示提示:非线性最小二乘问题非线性最小二乘问题,最速下降法或最速下降法或Gauss-Newton方法求解方法求解40/4614:50例例4Google街景技

27、术关键部分街景技术关键部分(大型复杂的非线性最小大型复杂的非线性最小二乘问题二乘问题)http:/google-opensource.blogspot.ca/2012/05/introducing-ceres-solver-nonlinear.html41/4614:50最小二乘拟合问题研究包括最小二乘拟合问题研究包括:模型的选取模型的选取(解释为什么选这个模型解释为什么选这个模型)存在唯一性存在唯一性最小二乘解的计算最小二乘解的计算(线性与非线性线性与非线性)42/4614:50不同变形不同变形2221,argminniairAabrr 22argminaAay 210argmin()sq(

28、uares)mnikkiaikyaxleast Aay 的的最最小小二二乘乘解解010110()()()()nmmnmnayxxryxxa0,()nikkikyax 观观测测值值模模型型预预测测值值43/4614:50数学概念数学概念1秩亏缺秩亏缺(rank deficient)矩阵矩阵A属于属于Rm*n如果如果rank(A)=m(n),则称矩阵行满秩则称矩阵行满秩如果如果rank(A)=n(m),则称矩阵列满秩则称矩阵列满秩如果如果rank(A)minm,n,则称矩阵秩亏缺则称矩阵秩亏缺44/4614:50数学概念数学概念2如果函数如果函数f 的的f,f,.,f(k)存在且连续。存在且连续。

29、理论分析中总是假设观测数据是由给定的理论分析中总是假设观测数据是由给定的函数生成函数生成(假设函数的数学性态假设函数的数学性态),然后定量的然后定量的刻画插值函数或者拟合函数的近似程度刻画插值函数或者拟合函数的近似程度。Ck 函数类函数类45/4614:50数学概念数学概念3极值极值(the maximum and minimum,known collectively as extrema)within a given neighborhood(local extremum)or on the function domain entirely (global extremum)极值点极值点(maximum point and minimum point)46/4614:50数学概念数学概念3 2 argmin()argmin()=argmin(),0()0,argmin()argmin()xxxxxf xf xccf x cf xf xf x 进进一一步步如如果果无约束优化问题无约束优化问题中中某点为极值点的某点为极值点的必要条必要条件件是该点梯度等于零。是该点梯度等于零。进一步地如果目标函数为凸函数进一步地如果目标函数为凸函数,某点为极某点为极值点的值点的充分必要条件充分必要条件为该点梯度等于零。为该点梯度等于零。