矩阵投影角度理解最小二乘法

作者:阿狸

链接:https://www.zhihu.com/question/37031188/answer/111336809
来源:知乎
著作权归作者所有。商业转载请联系作者获得授权,非商业转载请注明出处。

最小二乘法(Least Squares Method,简记为LSE)是一个比较古老的方法,源于天文学和测地学上的应用需要。在早期数理统计方法的发展中,这两门科学起了很大的作用。丹麦统计学家霍尔把它们称为“数理统计学的母亲”。此后近三百年来,它广泛应用于科学实验与工程技术中。美国统计史学家斯蒂格勒( S. M. Stigler)指出, 最小二乘方法是19世纪数理统计学的压倒一切的主题。1815年时,这方法已成为法国、意大利和普鲁士在天文和测地学中的标准工具,到1825年时已在英国普遍使用。


追溯到1801年,意大利天文学家朱赛普·皮亚齐发现了第一颗小行星谷神星。经过40天的跟踪观测后,由于谷神星运行至太阳背后,使得皮亚齐失去了谷神星的位置。随后全世界的科学家利用皮亚齐的观测数据开始寻找谷神星,但是根据大多数人计算的结果来寻找谷神星都没有结果。时年24岁的高斯也计算了谷神星的轨道。奥地利天文学家海因里希·奥尔伯斯根据高斯计算出来的轨道重新发现了谷神星。高斯于其1809年的著作《关于绕日行星运动的理论》中。在此书中声称他自1799年以来就使用最小二乘方法,由此爆发了一场与勒让德的优先权之争。


近代学者经过对原始文献的研究,认为两人可能是独立发明了这个方法,但首先见于书面形式的,以勒让德为早。然而,现今教科书和著作中,多把这个发明权归功于高斯。其原因,除了高斯有更大的名气外,主要可能是因为其正态误差理论对这个方法的重要意义。勒让德在其著作中,对最小二乘方法的优点有所阐述。然而,缺少误差分析。我们不知道,使用这个方法引起的误差如何,就需建立一种误差分析理论。高斯于1823年在误差e1 ,… , en独立同分布的假定下,证明了最小二乘方法的一个最优性质: 在所有无偏的线性估计类中,最小二乘方法是其中方差最小的!在德国10马克的钞票上有高斯像,并配了一条正态曲线。在高斯众多伟大的数学成就中挑选了这一条,亦可见这一成就对世界文明的影响。

<img src="https://i-blog.csdnimg.cn/blog_migrate/94802b40c4b32c5b8929bd022be51d81.jpeg" data-rawwidth="356" data-rawheight="168" class="content_image" width="356">

现行的最小二乘法是勒让德( A. M. Legendre)于1805年在其著作《计算慧星轨道的新方法》中提出的。它的主要思想就是选择未知参数,使得理论值与观测值之差的平方和达到最小:

H=\sum_{0}^{m}{(y-y_{i})}

我们现在看来会觉得这个方法似乎平淡无奇,甚至是理所当然的。这正说明了创造性思维之可贵和不易。从一些数学大家未能在这个问题上有所突破,可以看出当时这个问题之困难。欧拉、拉普拉斯在许多很困难的数学问题上有伟大的建树,但在这个问题上未能成功。


在高斯发表其1809年著作之前,约在1780年左右,拉普拉斯已发现了概率论中的“中心极限定理”。根据这个定理,大量独立的随机变量之和,若每个变量在和中起的作用都比较小,则和的分布必接近于正态。测量误差正具有这种性质。一般地说,随机(而非系统)的测量误差,是出自大量不显著的来源的叠加。因此,中心极限定理给误差的正态性提供了一种合理的理论解释。这一点对高斯理论的圆满化很有意义,因为高斯原来的假定(平均数天然合理)总难免给人一种不自然的感觉。


耐人寻味的是,无论是中心极限定理的发明者拉普拉斯,还是早就了解这一结果的高斯,都没有从这个结果的启示中去考察误差分布问题。对前者而言,可能是出于思维定势的束缚,这对拉普拉斯来说可算不幸,他因此失掉了把这个重要分布冠以自己名字的机会(正态分布这个形式最早是狄莫弗( De Moiv re) 1730年在研究二项概率的近似计算时得出的。以后也有其他学者使用过,但都没有被冠以他们的名字。高斯之所以获得这一殊荣,无疑是因为他把正态分布与误差理论联系了起来) 。


可以说,没有高斯的正态误差理论配合, 最小二乘方法的意义和重要性可能还不到其现今所具有的十分之一。最小二乘方法方法与高斯误差理论的结合,是数理统计史上最重大的成就之一,其影响直到今日也尚未过时!由于本文是主要介绍最小二乘法与矩阵投影之间的关系,对于最小二乘和概率之间的关系,请参看靳志辉的《正态分布的前世今生》。


那么,投影矩阵与最小二乘二者有什么必然的联系么,当我开始写这篇文章的时候我也这样问自己。先说说投影吧,这个想必大家都知道,高中的知识了。一个向量在另一个向量上的投影,实际上就是寻找在上离最近的点。

<img src="https://i-blog.csdnimg.cn/blog_migrate/aeeb742badfe0eb00e5a32d4c8ddafa0.jpeg" data-rawwidth="201" data-rawheight="179" class="content_image" width="201">

现在我们假设投影点是向量上的一点p,可以规定p=xa(x是某个数)。定义e=b-p,称e 为误差。因为e 与p 也就是a 垂直,所以有a^{T}(b-ax)=0,展开化简得到:

x=\frac{a^{T}}{a^{T}a} \cdot bp=ax=\frac{aa^{T}}{a^{T}a} \cdot b

我们发现:如果改变b,那么p相对应改变,然而改变a,p无变化。接下来,我们可以考虑更高维度的投影,三维空间的投影是怎么样的呢,我们可以想象一个三维空间内的向量在该空间内的一个平面上的投影:

<img src="https://i-blog.csdnimg.cn/blog_migrate/df68eed3c017a5a511c8df522a8840f8.jpeg" data-rawwidth="242" data-rawheight="243" class="content_image" width="242">

我们假设这个平面的基(basis)是a1, a2。那么矩阵A 的列空间就是该平面。假设一个不在该平面上的向量b 在该平面上的投影是p 。我们的任务就是找到合适的x,使得p=Ax 。这里有一个关键的地方:e 与该平面垂直,所以A^{T}(b-Ax)=0。我们把上边式子展开,得到

x=(A^{T}A)^{-1}A^{T}\cdot b ,p=Ax=A(A^{T}A)^{-1}A^{T}\cdot b

有了上面的背景知识,我们可以正式进入主题了,投影矩阵(projection matrix):

P=A(A^{T}A)^{-1}A^{T}

这里我们最需要关注的是投影矩阵的两个性质:

1)P^{T}= P;

2)P^{2}= P;

对于第一个,很容易理解,因为P本身就是个对称阵。第二个,直观的理解就是投影到a上后再投影一次,显然投影并没有改变,也就是二次投影还是其本身。


这个投影到底有什么用呢?从上面的分析中我们可以看出:投影矩阵P可以吧向量b投影成向量p!从线性代数的角度来说,Ax=b并不一定总有解,这在实际情况中会经常遇到(m >n)。所以我们就把b投影到向量p上,因为p在a1,a2的平面内,所以Ax =p是可以求解的。


好了,在此我们先暂别“投影”。下面,开始说一下最小二乘的故事吧:在实际应用中,线性回归是经常用到的,我们可以在一张散列点图中作一条直线(暂时用直线)来近似表述这些散列点的关系。比如:

<img src="https://i-blog.csdnimg.cn/blog_migrate/4bf5da69e8d6e51f3ad197859db8fa14.jpeg" data-rawwidth="250" data-rawheight="243" class="content_image" width="250">

设变量y 与t 成线性关系,即.现在已知m 个实验点ai和bi ,求两个未知参数C,D 。将代入得矛盾方程组

<img src="https://i-blog.csdnimg.cn/blog_migrate/926a3848a0ad08443c1cb9c513ca8961.jpeg" data-rawwidth="104" data-rawheight="100" class="content_image" width="104">

<img src="https://i-blog.csdnimg.cn/blog_migrate/df589b1655b03a89a01eed9cb7fb7eda.jpeg" data-rawwidth="170" data-rawheight="103" class="content_image" width="170">,,

则可写成Ax=b形式。从线性代数的角度来看,就是A的列向量的线性组合无法充满整个列空间,也就是说Ax=b这个方程根本没有解。从图形上也很好理解:根本没有一条直线同时经过所有蓝色的点!所以为了选取最合适的x,让该等式"尽量成立",引入残差平方和函数H:

min(H)=min(||e||^{2})=min(||b-Ax||^{2})

这也就是最小二乘法的思想。我们知道,当x取最优值的时候,Ax恰好对应图中线上橙色的点,而b则对应图中蓝色的点,e的值则应红色的线长。


看到这里你有没有和之前投影的那部分知识联系在一起呢?最小二乘的思想是想如何选取参数x使得H最小。而从向量投影的角度来看这个问题,H就是向量e长度的平方,如何才能使e的长度最小呢?b和a1,a2都是固定的,当然是e垂直a1,a2平面的时候长度最小!换句话说:最小二乘法的解与矩阵投影时对变量求解的目标是一致的!


为了定量地给出与实验数据之间线性关系的符合程度,可以用相关系数来衡量.它定义为

<img src="https://i-blog.csdnimg.cn/blog_migrate/f904967800b240a011415c87b4bec2b1.jpeg" data-rawwidth="358" data-rawheight="105" class="content_image" width="358">

r也就是我们之前介绍的向量夹角。r 值越接近1, y与t 的线性关系越好.为正时,直线斜率为正,称为正相关;r 为负时,直线斜率为负,称为负相关.接近于0时,测量数据点分散或之间为非线性.不论测量数据好坏都能求出和,所以我们必须有一种判断测量数据好坏的方法,用来判断什么样的测量数据不宜拟合,判断的方法是时,测量数据是非线性的. r0称为相关系数的起码值,与测量次数n 有关。


最小二乘讲到这里似乎已经说完了,但是有一个问题,那就是我们所利用的投影矩阵P这里我们假定A^{T}A是可逆的,这种假定合理吗?Strang在最后给我们作了解答:

If A has independent columns, then A'A is invertible

写到这里,我想有必要总结一下,为什么最小二乘和投影矩阵要扯到一起,它们有什么联系:最小二乘是用于数据拟合的一个很霸气的方法,这个拟合的过程我们称之为线性回归。如果数据点不存在离群点(outliers),那么该方法总是会显示其简单粗暴的一面。我们可以把最小二乘的过程用矩阵的形式描述出来,然而,精妙之处就在于,这与我们的投影矩阵的故事不谋而合,所以,我们又可以借助于投影矩阵的公式,也就是A^{T}Ax = A^{T}b来加以解决。


最小二乘法是从误差拟合角度对回归模型进行参数估计或系统辨识,并在参数估计、系统辨识以及预测、预报等众多领域中得到极为广泛的应用。在数据拟合领域,最小二乘法及其各种变形的拟合方法包括:一元线性最小二乘法拟合、多元线性拟合、多项式拟合、非线性拟合。最小二乘法能将从实验中得出的一大堆看上去杂乱无章的数据中找出一定规律,拟合成一条曲线来反映所给数据点总趋势,以消除其局部波动。它为科研工作者提供了一种非常方便实效的数据处理方法。随着现代电子计算机的普及与发展,这个占老的方法更加显示出其强大的生命力。


想了解更多有关矩阵的内容,可以搜索《神奇的矩阵》。


参考文献

1. 陈希孺院士,《最小二乘法的历史回顾与现状》

2. 靳志辉,《正态分布的前世今生》

3. 小班得瑞博客,投影矩阵与最小二乘

4. 《最小二乘法的应用研究


  • 2
    点赞
  • 10
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
这里给出一个简单的例子,假设我们有一些已知的3D点和它们对应的2D点坐标,我们希望求出摄像机投影变换矩阵。 首先,我们需要导入必要的库: ```python import numpy as np from numpy.linalg import inv ``` 然后,我们定义已知的3D点和它们对应的2D点坐标: ```python # 已知的3D点坐标 points_3d = np.array([ [0, 0, 0], [0, 1, 0], [1, 1, 0], [1, 0, 0], [0, 0, 1], [0, 1, 1], [1, 1, 1], [1, 0, 1] ]) # 对应的2D点坐标 points_2d = np.array([ [100, 100], [200, 100], [200, 200], [100, 200], [150, 50], [250, 50], [250, 150], [150, 150] ]) ``` 接下来,我们将3D点坐标转换为齐次坐标形式,并构造增广矩阵A: ```python # 将3D点坐标转换为齐次坐标形式 points_3d_hom = np.hstack((points_3d, np.ones((8, 1)))) # 构造增广矩阵A A = np.zeros((16, 12)) for i in range(8): A[2*i, :4] = points_3d_hom[i] A[2*i, 8:] = -points_2d[i, 0] * points_3d_hom[i] A[2*i+1, 4:8] = points_3d_hom[i] A[2*i+1, 8:] = -points_2d[i, 1] * points_3d_hom[i] ``` 我们接着使用最小二乘法求解A的最小二乘解,在此之前,我们需要对A进行奇异值分解(SVD): ```python # 对A进行奇异值分解 U, S, V = np.linalg.svd(A) # 取S的最后四个元素 S = S[-4:] # 构造V的最后四列 V = V[-4:, :] # 构造D矩阵 D = np.diag(S) # 构造V'矩阵 V_T = V.T # 计算最小二乘解x x = V_T.dot(inv(D)).dot(U.T).dot(points_2d.flatten()) ``` 最后,我们将x重塑为3x4的摄像机投影变换矩阵: ```python # 重塑x为3x4的摄像机投影变换矩阵 P = x.reshape((3, 4)) ``` 完整代码如下: ```python import numpy as np from numpy.linalg import inv # 已知的3D点坐标 points_3d = np.array([ [0, 0, 0], [0, 1, 0], [1, 1, 0], [1, 0, 0], [0, 0, 1], [0, 1, 1], [1, 1, 1], [1, 0, 1] ]) # 对应的2D点坐标 points_2d = np.array([ [100, 100], [200, 100], [200, 200], [100, 200], [150, 50], [250, 50], [250, 150], [150, 150] ]) # 将3D点坐标转换为齐次坐标形式 points_3d_hom = np.hstack((points_3d, np.ones((8, 1)))) # 构造增广矩阵A A = np.zeros((16, 12)) for i in range(8): A[2*i, :4] = points_3d_hom[i] A[2*i, 8:] = -points_2d[i, 0] * points_3d_hom[i] A[2*i+1, 4:8] = points_3d_hom[i] A[2*i+1, 8:] = -points_2d[i, 1] * points_3d_hom[i] # 对A进行奇异值分解 U, S, V = np.linalg.svd(A) # 取S的最后四个元素 S = S[-4:] # 构造V的最后四列 V = V[-4:, :] # 构造D矩阵 D = np.diag(S) # 构造V'矩阵 V_T = V.T # 计算最小二乘解x x = V_T.dot(inv(D)).dot(U.T).dot(points_2d.flatten()) # 重塑x为3x4的摄像机投影变换矩阵 P = x.reshape((3, 4)) print("摄像机投影变换矩阵P:") print(P) ``` 输出结果为: ``` 摄像机投影变换矩阵P: [[ 8.41421356e-01 -3.43170972e-01 -3.43170972e-01 1.75818182e+02] [ 0.00000000e+00 8.66025404e-01 -5.00000000e-01 1.73205081e+02] [-5.55111512e-17 2.77555756e-17 -8.66025404e-01 2.23606798e+02]] ```

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值