10

最小二乘法的参数估计

 3 years ago
source link: http://lanbing510.info/2016/03/28/Least-Squares-Parameter-Estimation.html
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.
neoserver,ios ssh client



计算机视觉中很多问题都涉及到模型参数估计,例如给定两幅图像中的匹配点对来求取单应变换矩阵。最小二乘本质上是一种通过最小化均方误差来求取最优的模型参数的数据拟合方法。下文将依次介绍线性最小二乘、非线性最小二乘及它们的求解方法。

线性最小二乘法


假设有p个关于q个未知数的线性方程组:

image01.jpg

在这个方程组中令

image02.jpg

由线性代数理论可知:

  • 当p<q时,方程的解集构成一个(q−p)维数的Rq的子空间;

  • 当p=q时,有唯一解;

  • 当p>q时,方程无解。

当U的秩(线性独立的行或列的最大数目)最大(等于min(p,q))时,以上结论成立。当秩小于min(p,q)时,解的存在性取决于y的取值以及它是否在U构成的空间内(由各列张成的Rp的子空间)。

正则方程和伪逆

本节剩余部分将假设U是过约束的,即p>q,且秩为q。在这种情况下没有标准的解,只能找到使误差最小的向量x,误差定义为

image03.jpg

上式最小化均方误差即为最小二乘法。将误差写为E=e⋅e,其中edef=Ux−y。为找到使E最小的x,xd 每一个分量的导数必须为0,即

image04.jpg

若U的各列表示为cj=(u1j,⋯,umj)T(j=1,⋯,q),则有

image05.jpg

带入前式有

image06.jpg

当U满秩时,矩阵UTU是可逆的,则方程x=U†y(其中U†def=(UTU)−1UT)的解就是所求的x。U†称为伪逆矩阵。在求解线性最小二乘法时,不一定要计算伪逆矩阵的实际值,一般通过QR或者SVD分解的方法来更快速的求解。具体求解步骤可以参考文献[2]。

齐次系统的求解

如果y为0,则原始问题转换为了齐次系统的求解问题,即

image07.jpg

上式是关于x的齐次方程。当p=q且矩阵U非奇异时,上式有唯一解;当p≥q且U奇异时(秩小于q),方程才有非零解。此种情况下如果不引入附加限制,原始的误差最小化E=|Ux|2是没有意义的,因为总能让x接近0来使得误差很小。由方程的齐次性,有E(λx)=λ2E(x),因此要求|x|2=1是一个合理限制,避免没有意义的解并且可以得到唯一解。

误差E可以写成|Ux|2=xT(UTU)x,其中q×q的矩阵UTU是半正定的,可以通过特征值分解过程对角化,特征向量为ei(i=1,⋯,q),特征值为0≤λ1≤⋯≤λq。任意单位向量x可用这些特征向量表示为x=u1e1+⋯+uqeq,其中u21+⋯+u2q=1,则有:

image08.jpg

可以看出使E最小的x即为UTU的最小特征值对应的特征向量e1,此时得到最小误差λ1。计算对称矩阵依旧可以用QR分解,SVD分解等多种方法,而不需要算出UTU。

非线性最小二乘法


考虑一般情况下的p个方程和q个未知数:

image09.jpg

其中fi(i=1,⋯,p)是Rp到R上的可微函数。用缩写f=(f1,⋯,fp)T和x=(x1,⋯,xq)表示。一般情况下有

  • 当p<q时,解构成一个Rq上的(q−p)维的子集;

  • 当p=q时,有有限个解;

  • 当p>q时,方程无解。

这与线性情况有几处明显的不同:在欠约束的情况下,解的维数仍是q−p,但是不再构成一个向量空间,其结构由fi决定。在p=q的情况下,不在是唯一解而是若干个。没有一个通用的方法可以找到上述方程在p=q的所有解或者p>q是的全局最小值解

image11.jpg

因此下面要介绍的是把原问题简化为局部线性问题的迭代方法,来力求找到至少一个合理解。它们都是建立在x领域上fi的一阶泰勒展开式基础上:

image12.jpg

其中▽fi(x)=(∂fi/∂x1,⋯,∂fi/∂xq)T是fi在点x的梯度值。在忽略二次项(O(|δx|2))的情况下马上可以得到

image13.jpg

其中Jf(x)是f的Jacobian矩阵,定义为一个p×q矩阵

image14.jpg

牛顿法:非线性方程的方阵系统

思想是:当p=q是,一般方法不能找到所有的解,但可以以f(x+δx)≈f(x)+Jf(x)δx为基础,构造一个迭代方法找到其中一个解。已知解的当前估计值为x,对x施加一个扰动δx来使得f(x+δx)=0,即

image15.jpg

当这个Jacobian矩阵非奇异时,解上面的方程既可以找到δx的合适解,重复这个过程直至收敛。牛顿法在接近解的地方收敛很快,按照平方速度收敛,但起始点离解很远时,上述牛顿法可能失败。

牛顿法:非线性方程的过约束系统

思想是:当p>q时,找到一个均方误差E的局部极小解。在极小值点,E的导数为零,利用这个特征来使用牛顿法。令F(x)=12▽E(x),用牛顿法可以找到非线性方程组F(x)=0的一组解。由E微分可得到

image17.jpg

F的Jacobian矩阵就是:

image20.jpg

其中Hfi(x)是fi的Hessian矩阵,由fi的二阶导数组成

image19.jpg

在牛顿法中,δx满足JF(x)δx=−F(x),则可以推出

image21.jpg

高斯牛顿法和Levenberg-Marquardt方法

牛顿法需要计算Hessians矩阵,不但困难而且耗时。下面介绍两种不需要计算Hessians矩阵的非线性优化方法:高斯牛顿法和Levenberg-Marquardt方法。

高斯牛顿法还是利用f的一阶泰勒展开式逼近E的极小值。但对给定的x,找到δx使得E(x+δx)最小。

image22.jpg

现在问题转化为了求线性最小二乘解的问题,δx可以通过解线性方程或通过伪逆得到

image24.jpg

比较上式和非线性方程过约束系统的牛顿法可以看出高斯牛顿法是牛顿法的一种近似,其忽略了Hessians矩阵部分。当fi在解的附近取值(残差)很小是这种近似是可以的。若残差很大,高斯牛顿法可能收敛很慢或者不收敛。

对上式稍加修改有

image25.jpg

其中μ在每步迭代可以取不同的值,此即为计算机视觉中常用的Levenberg-Marquardt方法。其用单位阵的倍数取代了包含Hessians的项。其和高斯牛顿法有相同的收敛速度,但是它更加鲁棒:即时在Jacobian矩阵Jf不满秩以及伪逆不存在的情况下也能够使用。


上文基本摘录整理自《计算机视觉-一种现代方法》第三章,感兴趣的同学根据参考文献进行延伸阅读。


[1] Computer Vision: A Modern Approach. David A. Forsyth and Jean Poince.

[2] 求解最小二乘的几种方法

[3] 计算机视觉-一种现代方法. David A. Forsyth and Jean Poince.


About Joyk


Aggregate valuable and interesting links.
Joyk means Joy of geeK