数值分析16最小二乘法

1766年，德国有一位名叫约翰·提丢斯的中学教师发现如下数列的每一项与当时已知的六大行星(即水星、金星、地球、火星、木星、土星)到太阳的距离比例(地球到太阳的距离定为1个单位)有着一定的联系(0+4)/10=0.4；(3+4)/10=0.7；(6+4)/10=1.0；(12+4)/10=1.6；

(24+4)/10=2.8；(48+4)/10=5.2；

(96+4)/10=10；(3*2ⁿ+4)/10？1781年，英籍德国人赫歇尔在接近19.6的位置上(即数列中的第八项)发现了天王星，从此，人们就对这一定则深信不疑了。根据这一定则，在数列的第五项即2.8的位置上也应该对应一颗行星，只是还没有被发现。于是，许多天文学家和天文爱好者便以极大的热情，踏上了寻找这颗新行星的征程。1801年新年的晚上，意大利神父皮亚齐还在聚精会神地观察着星空。突然，他从望远镜里发现了一颗非常小的星星，正好在提丢斯——波得定则中2.8的位置上。这颗行星在几天的观测期内不断变动位置。可是，当皮亚齐再想进一步观察这颗小行星时，他却病倒了。等到他恢复健康，再想寻找这颗小行星时，它却不知所踪。高斯根据皮亚齐的观测资料，利用最小二乘方法，只用了一个小时就算出了谷神星的轨道形状，并指出它将于何时出现在哪一片天空里。1801年12月31日夜，德国天文爱好者奥伯斯，在高斯预言的时间里，用望远镜对准了这片天空。不出所料，谷神星再一次奇迹般地出现了！高谷神星谷神星斯的计算方法成功了。高斯从笔尖上寻找到的这颗行星，在隐藏了整整一年后，成为人类的最好的新年礼物。这一礼物向人们显示了数学在科学研究中的巨大作用。无解的方程组称为不相容方程组。方程个数多于未知数个数的方程组称为超定方程组。定义:一个方程组若至少存在一个解能够严格满足该方程组,则称为相容方程(consistentequation)。定理:线性方程Ax=b相容的充分必要rank([A,b])=rank(A)

Thehardestquantitiestocomputeareonesthatdonotexist.

当方程组不相容时,如何寻求次佳(next-best)解。我们希望残差r=Ax-b尽可能小,小的标准是什么呢？最小二乘解

(least-squares)别称：拟合(fitting)、回归(regression)多变量函数极值点:几何解释？满足插值条件(1,2),(-1,1)和(1,3)的一次多项式。最小二乘解：正规方程：不相容方程的最小二乘解总是存在的。证明:即证明正规方程是相容方程。

rank([ATA,b])=rank(ATA)rank([ATA,b])=rank(ATA)相容

不相容唯一解不唯一解(最小范数解)最小二乘解参考:张贤达,矩阵分析与应用,

清华大学唯一解不唯一解(最小范数二乘解)例1:美国奥勒冈州波特兰市47套房产的售价和居住面积数据如下:Ref：求拟合(fitting)函数:

x

x1

x2··········xmyy1

y2··········ym离散数据的拟合求拟合函数:正规方程：

x

x1

x2··········xmy

y1

y2··········ym离散数据的拟合求拟合函数:正规方程：最小二乘拟合的步骤:给定的离散观测数据(x1,y1),…,(xm,ym)1.选取模型:确定模型类型(或拟合函数),如y=a+bx。2.将数据代入模型,每个观测数据对应一个参数作为未知数的方程。3.求解最小二乘解(即确定模型的参数)。观测数据(也许数据背后没有规律)发现规律解释规律背后的机理报道称，胡立德教授与其团队使用一架高速摄影机拍摄包括老鼠、山羊、牛、大象、狗等哺乳动物的排尿过程，并通过流体动力学建模发现，只要体重超过3公斤的哺乳动物，它们的排尿时间都在21秒左右。例.从20世纪70年代早期开始,英特尔公司CPU上的晶体管的数量:CPUyeartransistors400419712,250800819722,500808019745,0008086197829,0002861982120,0003861985275,00048619891,180,000Pentium19933,100,000Pentium

II19977,500,000Pentium

III199924,000,000Pentium

4200042,000,000Itanium2002220,000,000Itanium

22003410,000,000指数函数刻画增长或衰减的规律,如放射性原子核、血液中药物和酒精的衰变规律和。摩尔定律是由英特尔（Intel）创始人之一戈登·摩尔（GordonMoore）提出来的。其内容为：当价格不变时，集成电路上可容纳的晶体管数目，约每隔18(或24)个月便会增加一倍，性能也将提升一倍。线性化改变了原最小二乘问题。对于原非线性最小二乘问题可以考虑最速下降法和共轭梯度法等优化方法。网络科学领军人物艾伯特-拉斯洛·巴拉巴西写过的两本书《链接》和《爆发》。他在这两本书里总结出的是，世界是由两个法则构成的，一个是高斯法则，也就是正态分布，简单说就是“中间多，两头少”，比如一个班里，成绩中间的学生最多，成绩特别好和成绩特别差的学生不会太多。另一个是幂律法则，可以理解为二八法则，比如20%的人口占据了80%的社会财富。幂律行星规律探秘行星轨道半长轴(1010米)周期(年)水星5.790.241金星10.80.615地球15.01火星22.81.88木星77.811.9土星14329.5天王星28784海王星450165冥王星590248

开普勒大胆猜测：“绕以太阳为焦点的椭圆轨道运行的所有行星,其各自椭圆轨道半长轴的立方与周期的平方之比是一个常量?行星轨道半长轴(1010米)周期(年)水星5.790.241金星10.80.615地球15.01火星22.81.88木星77.811.9土星14329.46天王星28784.01海王星450164.8冥王星590247.7线性化改变了原最小二乘问题,对于原非线性最小二乘问题可以考虑最速下降法和

共轭梯度法等优化方法。

生物界的牛顿定律物種面積曲線幂律MaxKleiber在1930年发现了一条在整个生物界都普遍适用的定律：一个动物的能耗n与它体重的关系为n=v^(3/4)。如何选择最佳的模型?若无必要，勿增实体奥卡姆剃刀原理:当两个假说具有完全相同的解释力和预测力时，我们以那个较为简单的假说作为讨论依据。

主要函数1.多项式曲线拟合函数：polyfit()调用格式： p=polyfit(x,y,n)

说明：x,y为数据点，n为多项式阶数，返回p为幂次从高到低的多项式系数向量p。2.多项式曲线求值函数：polyval()调用格式： y=polyval(p,x)说明：y=polyval(p,x)为返回对应自变量x在给定系数P的多项式的值。x=[1245];y=[1223];plot(x,y,'o')p=polyfit(x,y,1)xnew=linspace(min(x),max(x));ynew=polyval(p,xnew);holdon,plot(xnew,ynew,'r')研究美国人口增长的规律并预测2020年的人口。Year19001910192019301940

1950Population75.99591.972105.711123.203131.669150.697Year196019701980199020002010Population179.323

203.212

226.505249.633281.422308.746loadcensus,plot(x,y,'o')P=polyfit(x,y,3);polyval(P,2010)Curve

FittingTool严格满足插值条件vs

追求最小残差平方和适定方程组求解vs

不相容方程组求解寻找数据的规律(或压缩数据)思考:插值与(最小二乘)拟合的异同？vs

预测的能力作业5研究饮酒后人体血液中酒精含量的变化规律Ref:时间0.250.50.7511.522.533.544.55酒精含量306875827768686858515041时间678910111213141516酒精含量3835282518151210774体重约70kg的某人在短时间内喝下2瓶啤酒后，隔一定时间测量他的血液中酒精含量(毫克／百毫升)，数据如下：作业6Sunspots

Ref:C.Moler,NumericalComputingwithMatlabForcenturies,peoplehavenotedthatthefaceofthesunisnotconstantoruniforminappearance,butthatdarkerregionsappearatrandomlocationsonacyclicalbasis.Thisactivityiscorrelatedwithweatherandothereconomicallysignificantterrestrialphenomena.In1848,RudolfWolferproposedarulethatcombinedthenumberandsizeofthesesunspotsintoasingleindex.Usingarchivalrecords,astronomershaveappliedWolfer’sruletodeterminesunspotactivitybacktotheyear1700.

微信:狗熊会(北京大学王汉生教