码迷,mamicode.com
首页 > 其他好文 > 详细

漫步线性代数十六——投影和最小二乘

时间:2016-09-08 18:35:57      阅读:245      评论:0      收藏:0      [点我收藏+]

标签:

目前为止,我们已经知道Ax=b要么有解要么无解,如果b 不在列空间C(A) 里,那么这个系统就是矛盾的,高斯消元法就会失败。当有几个方程和一个未知量时失败完全可以确定:

2x3x4x===b1b2b3

b1,b2,b3的比率是2:3:4时,上面的方程组才可解,也就是说只有b 和列a=(2,3,4)在一条直线上时x才会存在。

尽管他们无解,可是他们在实际中经常出现,他们必须有解!一种可能是用系统的一部分来确定x,其余部分忽略;如果所有的m个方程来源一样,这种方法就不合理。我们放弃这种一些方程没误差,而有些误差大的想法,我们考虑能最小化m个方程平均误差Ex值。

对平方和求平均是最方便的:

E2=(2x?b1)2+(3x?b2)2+(4x?b3)2

如果存在准确解,那么最小误差E=0。大部分情况下,ba不成比例关系,E2的图像将是一个抛物线,最小误差在最低点的位置处,也就是导数等于零的位置:

dE2dx=2[(2x?b1)2+(3x?b2)3+(4x?b3)4]=0

求出x的值,这个模型系统ax=b的最小二乘解用x^ 来表示:

x^=2b1+3b3+4b322+32+42=aTbaTa

相信大家立马就认出分子中的aTb和分母中的aTa了吧(是不是像投影啊)。

推广到一般情况同样如此,求解ax=b就是最小化

E2=ax?b2=(a1x?b1)2+?+(amx?bm)2

E2求导并令其等于零,求出点x^

(a1x^?b1)a1+?+(amx^?bm)am=0

计算后得到x^=(a1b1+?+ambm)/(a21+?+a2m)

11、对于ax=b这样只有一个未知变量的问题,它的最小二乘解为:x^=aTbaTa

大家可能看出来了,我们一直从几何角度解释最小二乘问题—— 最小化距离。令E2的导数等于零求出解,求得的结果和上篇文章的几何形式一样,连接b,p 的误差向量e一定垂直于a

aT(b?x^a)=aTb?aTbaTaaTa=0

注意退化为a=0的情况,这是a的任何倍数都是零,线仅仅就是一个点,因此p=0是唯一的投影候选结果。但是x^的形式变成一个无意义的数0/0,这表明x^完全无法确定,所有x值都给出相同的误差E=0x?b,所以E2是一条水平线而不是抛物线,伪逆给这种情况分配了一个确定的值x^=0,相比较其他值,这个是最好的选择的。

最小二乘问题

现在我们开始难一点的问题,将b投影到一个子空间上——而不是一条线上。这个问题来自于Ax=b,其中Am×n矩阵,不再是一列和一个未知量x,现在矩阵有多列,m 的个数比未知量n的个数要大,所以跟期望中的一样,Ax=b依然是矛盾的。不可能存在完全拟合数据bx值,换句话说,向量b不是A列向量的组合;它在列空间的外面。

再次回到了找出x^来最小化误差的问题,这个最小化可以用最小二乘求解,误差是E=Ax?b,这就是b到列空间中Ax的距离。我们要做的就是能最小化E的最小二乘解x^,它和p=Ax^相等,而这个p就是列空间中离b最近的点。

我们可以用几何或计算来确定x^,在n维空间中,我们偏爱几何;p 一定是b在列空间上的投影。误差向量e=b?Ax^一定可这个空间垂直(图1),找到x^和投影p=Ax^是最基本的,下面我们用两种方法来实现它:

  1. 所有垂直于列空间的向量位于左零空间里,因此误差向量e=b?Ax^一定在AT的零空间里:
    AT(b?Ax^)=0orATAx^=ATb
  2. 误差向量和A的每列a1,,an垂直:
    aT1(b?Ax^)=0?aTn(b?Ax^)=0or????aT1?aTn????[b?Ax^]=0


    技术分享
    图1

这两种方法殊途同归,最后都是AT(b?Ax^)=0,ATAx^=ATb,而计算方法是通过计算E2=(Ax?b)T(Ax?b)的导数,并令其等于零得2ATAx?2ATb=0,最快的方式是方程Ax=b两边乘以AT,所有这些等价方法都得到一个二次系数矩阵ATA,它是对称的(它的转置可不是AAT!)并且是接下来几篇文章中非常基础的矩阵。

方程ATAx^=ATb在统计学中叫做正规方程。

12、当Ax=b是矛盾的时候,它的最小二乘解就是最小化Ax?b2

ATAx^=ATb(1)
A的列线性无关时,ATA是可逆的!因此
x^=(ATA)?1ATb(2)
b在列空间上的投影就是最近点Ax^
p=Ax^=A(ATA)?1ATb(3)

我们举一个例子进行说明:

A=???110230???,b=???456???,Ax=b,ATAx^=ATbx

每个列最后一个元素都是零,所以C(A)是三维空间中的x?y平面,b=(4,5,6)的投影是p=(4,5,0)x,y分量保持不变,但z分量变成零,通过求解正规方程就能证实这个结果:

ATA=[121300]???110230???=[25513]

x^=(ATA)?1ATb=[13?5?52][121300]???456???=[21]

p=Ax^=???110230???[21]=???450???

在这种特殊情况,最佳方式就是求解Ax=b的前两个方程,得到x^1=1,x^2=1,方程0x1+0x2=6的误差是6。

注解:假设bA的列空间里,也就说存在列的组合使得b=Ax,那么b的投影依然是b

p=A(ATA)?1ATAx=Ax=b

最近的点p就是b本身。

注解:考虑一个极端的情况,假设b与每列都垂直,那么ATb=0,这种情况下b的投影就是零向量:

p=A(ATA)?1ATAx=A(ATA)?10=0

注解:当A是方阵且可逆时,列空间就是整个空间,每个向量的投影就是自身,p=b,x^=x

p=A(ATA)?1ATAx=AA?1(AT)?1ATb=b

只有这一种情况我们可以将(ATA)?1分离成A?1(AT)?1,当A是长方形矩阵时,就不能这么做。

注解:假设A只有一列,也就是只包含a,那么矩阵ATA就是常数aTax^就是aTb/aTa,回到了最初的形式。

矩阵ATA

矩阵ATA一定是对称的,因为它的转置(ATA)T=ATATT,依然是ATA。它的第i,j(j,i) 个元素是A的第i列与第j行的内积,重点是ATA 的可逆性,幸运的是ATAA有相同的零空间。如果Ax=0,那么ATAx=0A零空间中的向量x也在ATA的零空间中。反过来考虑,假设ATAx=0,我们将它和x进行内积操作来表明Ax=0

xTATAx=0,orAx2=0,orAx=0

两个零空间是相等的。如果A有无关列(零空间中只有x=0),那么ATA同样如此:

13、如果A有无关列,那么ATA是方阵,对称并且可逆。

随后我们还会指出ATA也是正定的(所有主元和特征值都是正的)。

到目前为止,这种情况是最常见也是最终要的,如果m>n,那么m维空间的无关性就很容易实现。

投影矩阵

我们已经说明了离b的最近点是p=A(ATA)?1ATb,这种形式用矩阵形式来表示就是构建bA列空间的垂线,产生p的矩阵是一个投影矩阵,用P 表示:

P=A(ATA)?1AT(4)

这个矩阵将任何向量b投影到A的列空间上,换句话说,p=Pbb在列空间上的分量,误差e=b?Pb是正交补中的分量。(I?P也是一个投影矩阵!它将b投影到正交补上,投影是b?Pb)

简单来说,有一种矩阵形式可以将b分成两个互相垂直的分量,Pb在列空间C(A)内,其他的分量(I?P)b在左零空间N(AT)内——也就是与列空间正交的空间。

这些投影矩阵可以从代数和几何两个角度理解。

14、投影矩阵P=A(ATA)?1AT有两个性质:

- 矩阵等于自身的平方:P2=P
- 矩阵等于它的转置:PT=P

反过来讲,任何对称矩阵,如果P2=P,那么它表示一种投影。

证明:很容易看出来为什么P2=P,我们先从任意向量b开始,那么Pb位于投影的子空间内,当我们再次投影的话不会发生任何变化,向量Pb已经在子空间内,P(Pb)依然是Pb,换句话说P2=P,两次或三次或五次投影得到的结果跟第一次一样:

P2=A(ATA)?1ATA(ATA)?1AT=A(ATA)?1AT=P

为了证明P是对称的,我们取它的转置:

PT=(AT)T((ATA)?1)TAT=A(ATA)?1AT=P

反过来,我们可以从P2=P,PT=P推断出PbbP列空间上的投影,误差向量b?Pb与这个空间正交。对于该空间内的所有向量Pc,内积是零:

(b?Pb)TPc=bT(I?P)TPc=bT(P?P2)c=0

因此b?Pb和空间是正交的,Pb是列空上的投影。

例1:假设A是可逆的,如果它是4×4矩阵,那么它的四列都是无关的,列空间就是整个R4。在整个空间上的投影是什么?答案就是单位矩阵。

P=A(ATA)?1AT=AA?1(AT)?1AT=I(5)

单位矩阵是对称的,并且I2=I,误差向量b?Ib等于零。

拟合数据的最小二乘法

假设我们有一堆实验数据,并且期望输出b是输入t的线性函数,也就是看成直线b=C+Dt,例如:

  1. 我们测量不同时刻卫星距火星的距离,我们用t表示时间,b表示时间,不考虑失去动力或重力突然增强的情况下,卫星几乎以恒定的速度v移动:b=b0+vt
  2. 我们在某个物体上放上不同的载荷,并测量它垂直方向产生的位移,我们用t 表示载荷的重量,b表示位移大小。除非载太重使得物体彻底变形,否则的话根据弹性理论,存在一个线性关系b=C+Dt
  3. 印制t本书的成本似乎也是线性关系:b=C+Dt,其中编辑和排版成本是C,印刷和装订成本是DC是固定的,而每印制一本书成本多D

如何计算C,D呢?如果没有实验误差,那么两次测量的b都会得到直线b=C+Dt,但是如果有误差的话,我们就考虑平均值,求出最佳的直线。事实上,因为有两个未知量C,D需要确定,于是我们需要投影到二维子空间上。而一般情况下,我们都是多次进行试验测量的:

CCC+++Dt1Dt2?Dtm===b1b2bm(6)

得到的是矛盾方程组,有m个方程却只有两个未知量,如果误差存在的话,它将不可解。我们写成矩阵形式:

??????11?1t1t2?tm??????[CD]=??????b1b2?bm??????,orAx=b(7)

最佳解(C^,D^)就是最小化均方误差E2得到的x^

E2=b?Ax2=(b1?C?Dt1)2+?+(bm?C?Dtm)2

向量p=Ax^是最接近向量b的,在所有的直线b=C+Dt中,我们选出拟合数据最好的直线(图2),在图中,误差是到直线的竖直距离b?C?Dt(不是垂直距离!),它对应的是竖直距离的平方,求和和最小化。

例2:在图2a中有三个测量值b1,b2,b3

t=?1,b=1;t=1,b=1;t=2,b=3

注意t=?1,1,2不要求等距离。第一步是通过三个点的方程:

Ax=bisCCC?+?DD2D===113???111?112???[CD]=???113???

如果这些方程Ax=b可解,那么表示没有误差。但是这些点不在一条直线上,所以他们不可解,因此需要用到最小二乘求解:

ATAx^=ATb[3226][C^D^]=[56]

最佳解就是C^=97,D^=47,最佳直线是97+47t


技术分享
图2

注意这两幅图之间的联系,问题是一样的但是呈现的效果不一样。在图2b中,b不是列(1,1,1),(?1,1,2)的一个组合,而在图2a中,三个点不在一条线上。最小二乘用点p代替了不在直线上的点b!既然无法解Ax=b,那我们就解Ax^=p

直线97+47t?1,1,2处的高度分别为57,137,177,这些点都在之直线上,因此向量p=(57,137,177)在列空间里,而这个向量就是投影。图2b展示的是三维空间效果(如果有m个点就是m维)而图2a 是二维空间的效果(如果有n 个参数就是n维)。

b中减去p得到误差e=(27,?67,47),在图2a中就是竖直向量,他们是图2b中虚线向量的元素,这个误差向量与第一列(1,1,1)正交,因为?27+67+47=0,跟第二列也正交,所以它与列空间正交,属于左零空间。

问题:如果测量结果b=(27,?67,47)就是误差,那么最佳直线和解x^是什么呢?答案是:零,也就是水平轴,x=0^,投影是零。

我们总结一下拟合直线的方法,A的第一列包含1,第二列包含t,因此ATA包含1,t,t2的和:

15、给定点t1,?,tm处的测量值b1,?,bm,那么最小二乘求E2得到的直线C^+D^t为:

ATA[D^D^]=ATb[mΣtiΣtiΣt2i][C^D^]=[ΣbiΣtibi]

注解:最小二乘法不限于用直线拟合数据,在许多实验中关系不一定是线性的。假设我们有一些放射性材料,在不同时刻t可以通过仪器读出放射量b。现在我们知道这些材料是两种化学物质的混合物,还知道他们的半衰期(或衰减率),但是不知道每种的含量。如果我们用C,D 表示这两个未知量,那么仪器的结果更像是两个指数之和(不是直线):

b=Ce?λt+De?μt(8)

而实际测量中,仪器的结果存在误差,所以我们多测几次,分别在t1,,tm时刻测得b1,,bm,利用方程(8)近似满足:

Ax=bCe?λt1Ce?λtm++De?μt1?De?μtmb1bm

如果记录的次数超过两次m>2,那么我们可能无法求解,但是最小二乘原则将给出最佳解C^,D^

知道了C,D后情况就完全不同了,接下来我们就能算出衰减率λ,μ。这个问题就是非线性最小二乘,比线性的难一点。而我们依然是先写出E2,误差的平方和,然后最小化。但是导数为零得到的不再是线性方程。

加权最小二乘

一个简单的最小二乘问题是估计两个观测值x=b1,x=b2x^,除非b1=b2,否则我们面对的就是两个方程一个未知量的矛盾方程:

[11][x]=[b1b2]

目前为止,我们认为b1,b2可靠度一样,基于此我们最小化E2求出x^的值:

dE2dx=0x^=b1+b22

最佳解就是平均值,利用ATAx^=ATb得到同样的结果。事实上,ATA1×1的矩阵,正规方程是2x^=b1+b2

现在假设两个观测值的信任程度不一样,x=b1的结果比x=b2更加准确,但不管怎样,只要b2包含了信息,我们不会完全依赖b1,最简单的分解就是给他们分配不同的权值w21,w22,最下化带权重的平方和:

E2=w21(x?b1)2+w22(x?b2)2

如果w1>w2,那么说明b1更加重要,最小化过程时会使(x?b1)2变小的力度加大:

dE2dx=2[w21(x1?b1)+w22(x?b2)]=0,x^W=w21b1+w22b2w21+w22(9)

结果不再是b1,b2的平均值,而是数据的加权平均,这个平均相比b2更加靠近b1

一般最小二乘问题将Ax=b变成新系统WAx=Wb,这将结果x^变成了x^W,矩阵WTW出现在正规方程的两边:

WAx=Wb的最小二乘解是x^W

(ATWTWA)x^W=ATWTWb

b投影到Ax^的图像中发生了什么了?投影Ax^W依然是列空间中最靠近b的点,但是这里的最靠近有了新的意义,x的加权长度等于Wx的长度,垂直也不再是yTx=0,在新的方程组中是(Wy)T(Wx)=0,中间出现了矩阵WTW,在这个新观念下,投影Ax^W和误差b?Ax^W依然是垂直的。

接下里我们描述一下内积:他们来自于逆矩阵W。他们只涉及对称组合C=WTWx,y的内积是yTCx。对于正交矩阵W=Q,当这个组合是C=QTQ=I时,这和我们之前介绍的内积是一个含义,这种情况下旋转空间不改变内积,而其他矩阵会改变长度和内积。

对任何可逆矩阵W,这些规则定义了新的内积和长度:

(x,y)W=(Wy)T(Wx)xW=Wx(10)

因为W是可逆的,所以没有任何向量会变成零(除了零向量),所有可能的内积(线性依赖于x,y,并且在x=y0 时为正)可以从C=WTW 中找到。

实际中,重要的问题是C的选择,最好的答案来自统计学,最早是出自高斯。我们知道平均误差是零,这是b中误差的期望值(误差并非一定为零!),我们还知道误差平方的均值,也就是方差。如果bi的误差互相独立,且方差为σ2i,那么正确的权值是wi=1/σi,测量越精确(意味着更小的方差),权重越大。

除了不同的权重外,观测量也许是不独立的,如果误差是耦合的,那么W将是非对角形式,最好的非偏置矩阵C=WTW是协方差矩阵的逆(它的i,j项是bi误差和bj误差乘积的期望),C?1的主对角线包含方差σ2i,也就是bi误差平方的平均值。

例3:假设两个牌友(已经叫牌了)在猜对方手中黑桃的个数,误差为?1,0,1的概率都等于13,那么期望误差是零,方差是23

E(e)=13(?1)+13(0)+13(1)=0

E(e2)=13(?1)2+13(0)2+13(1)2=23

这两个人的猜测是相关的,因为叫牌是一样的,但是却不一样,这又是因为他们手中的牌不一样。如果说他们都猜大和都猜小的几率为零,相反误差的几率是13,那么E(e1e2)=13(?1),协方差矩阵的逆是WTW

[E(e21)E(e21e2)E(e21e2)E(e22)]?1=????23?13?1323?????1=[2112]?1=C=WTW

这就是加权正规方程中间的矩阵。

漫步线性代数十六——投影和最小二乘

标签:

原文地址:http://blog.csdn.net/u010182633/article/details/52444135

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!