RANSAC算法与原理（二）

文章目录

使用RANSAC进行直线拟合
总体最小二乘法直线推导
- 根据几何意义推导
- 从解齐次方程的角度推导
- - 方法一
  - 方法二
总结
References

之前介绍了RANSAC算法原理的基本原理和算法实现流程，那么我们在具体使用时，应当如何应用呢。由于RANSAC算法的通用性，其适用于各种模型拟合的场景下，如二维（三维）直线拟合、平面拟合、椭圆拟合等等。 skimage官方文档中给出了几个使用RANSAC进行模型拟合的例子。

使用RANSAC进行直线拟合

首先先看一个使用RANSAC进行直线拟合的简单例子。参考RANSAC算法详解(附Python拟合直线模型代码) - 知乎 (zhihu.com)，并对迭代次数的更新部分进行修改。

import numpy as np
import matplotlib.pyplot as plt
import random
import math# 数据量。
SIZE = 50
# 产生数据。np.linspace 返回一个一维数组，SIZE指定数组长度。
# 数组最小值是0，最大值是10。所有元素间隔相等。
X = np.linspace(0, 10, SIZE)
Y = 3 * X + 10fig = plt.figure()
# 画图区域分成1行1列。选择第一块区域。
ax1 = fig.add_subplot(1,1, 1)
# 标题
ax1.set_title("RANSAC")# 让散点图的数据更加随机并且添加一些噪声。
random_x = []
random_y = []
# 添加直线随机噪声
for i in range(SIZE):random_x.append(X[i] + random.uniform(-0.5, 0.5)) random_y.append(Y[i] + random.uniform(-0.5, 0.5))
# 添加随机噪声
for i in range(SIZE):random_x.append(random.uniform(0,10))random_y.append(random.uniform(10,40))
RANDOM_X = np.array(random_x) # 散点图的横轴。
RANDOM_Y = np.array(random_y) # 散点图的纵轴。# 画散点图。
ax1.scatter(RANDOM_X, RANDOM_Y)
# 横轴名称。
ax1.set_xlabel("x")
# 纵轴名称。
ax1.set_ylabel("y")# 使用RANSAC算法估算模型
# 迭代最大次数，每次得到更好的估计会优化iters的数值
iters = 100000
# 数据和模型之间可接受的差值
sigma = 0.25
# 最好模型的参数估计和内点数目
best_a = 0
best_b = 0
pretotal = 0
# 希望的得到正确模型的概率
P = 0.99
i=0
while True:if i<iters:i+=1# 随机在数据中红选出两个点去求解模型sample_index = random.sample(range(SIZE * 2),2)x_1 = RANDOM_X[sample_index[0]]x_2 = RANDOM_X[sample_index[1]]y_1 = RANDOM_Y[sample_index[0]]y_2 = RANDOM_Y[sample_index[1]]# y = ax + b 求解出a，ba = (y_2 - y_1) / (x_2 - x_1)b = y_1 - a * x_1# 算出内点数目total_inlier = 0for index in range(SIZE * 2):y_estimate = a * RANDOM_X[index] + bif abs(y_estimate - RANDOM_Y[index]) < sigma:total_inlier = total_inlier + 1# 判断当前的模型是否比之前估算的模型好if total_inlier > pretotal:iters = math.log(1 - P) / math.log(1 - pow(total_inlier / (SIZE * 2), 2))pretotal = total_inlierbest_a = abest_b = b# 判断是否当前模型已经符合超过一半的点if total_inlier > SIZE:breakelse:break
# 用我们得到的最佳估计画图
Y = best_a * RANDOM_X + best_b# 直线图
ax1.plot(RANDOM_X, Y)
text = "best_a = " + str(best_a) + "\nbest_b = " + str(best_b)
plt.text(5,10, text,fontdict={'size': 8, 'color': 'r'})
plt.show()

上述示例是使用RANSAC进行直线拟合的一个简单实现，我们可以注意到下述这段代码：

 for index in range(SIZE * 2):y_estimate = a * RANDOM_X[index] + bif abs(y_estimate - RANDOM_Y[index]) < sigma:total_inlier = total_inlier + 1

在进行内点和外点的判断选择时，是采用循环的方式逐一判断的，实际上，我们可以利用线性代数的知识，利用numpy矩阵运算将代码简化。

上述的代码，在每次进行直线拟合的时候，使用的两个点来进行直线方程的求解。那么当我想使用大于等于2个的数据点进行直线拟合时，最常见的就是使用最小二乘法拟合直线。

在使用最小二乘法进行直线拟合时，需要注意传统的最小二乘法（least squares)和总体最小二乘法(total least squares)的区别，详细说明请参考wiki链接，分别为最小二乘法 - 维基百科，自由的百科全书 (wikipedia.org)和[Total least squares - Wikipedia](https://en.wikipedia.org/wiki/Total_least_squares#:~:text=In applied statistics%2C total least squares is a,dependent and independent variables are taken into account.)。

两种方法的示意图如下图所示：

知乎上也有作者对这两个问题进行了说明和推导，参见线性拟合1-最小二乘 - 知乎 (zhihu.com)和

线性拟合2-正交回归 - 知乎 (zhihu.com)。

下面我将对总体最小二乘法的推导进行详细说明。

总体最小二乘法直线推导

根据几何意义推导

下面的推导主要是根据线性拟合2-正交回归 - 知乎 (zhihu.com)进行整理得到。

首先需要说明的是，下述推导是二维直线的点法式推导（三维直线没有点法式方程）。点法式的方程相比斜截式的方程，能更完整地的表示所有的直线情况。

横坐标和纵坐标的误差定义如下：
ηi^=xi−xi^\hat{\eta _{i}} =x_{i}-\hat{x_{i}} ηi^=xi−xi^

ϵi^=yi−yi^\hat{\epsilon _{i}} =y_{i}-\hat{y_{i}} ϵi^=yi−yi^

综合考虑横纵坐标的误差，目标函数的形式如下：
J2=∑in[(ηi^)2+(ϵi^)2]=∑in[(xi−xi^)2+(yi−yi^)2]J_{2}=\sum_{i}^{n} [(\hat{\eta _{i}})^{2}+(\hat{\epsilon _{i}})^{2}]=\sum_{i}^{n} [(x_{i}-\hat{x_{i}})^{2}+(y_{i}-\hat{y_{i}})^{2}] J2=i∑n[(ηi^)2+(ϵi^)2]=i∑n[(xi−xi^)2+(yi−yi^)2]
因为要求目标函数的最小值，所以就是求观测点(xi,yi)(x_{i},y_{i})(xi,yi)到估计点(xi^,yi^)(\hat{x_{i}},\hat{y_{i}})(xi^,yi^)的最小值，因为点(xi,yi)(x_{i},y_{i})(xi,yi)我们是已知的，也就是求点(xi,yi)(x_{i},y_{i})(xi,yi)到拟合直线的距离的和的最小值。估计点(xi^,yi^)(\hat{x_{i}},\hat{y_{i}})(xi^,yi^)是点(xi,yi)(x_{i},y_{i})(xi,yi)在拟合直线上的正交投影点。所以目标函数可以写成：
J2=∑indi2J_{2}=\sum_{i}^{n}d_{i}^{2} J2=i∑ndi2
did_{i}di为第iii个点(xi,yi)(x_{i},y_{i})(xi,yi)到拟合直线的距离。

用点法式拟合直线的表示形式为
a(x−x0)+b(y−y0)=0a(x-x_0)+b(y-y_0)=0 a(x−x0)+b(y−y0)=0
其中，(x0,y0)(x_0,y_0)(x0,y0)为经过直线的一个点，(a,b)(a,b)(a,b)表示直线的法向量。(a,b)(a,b)(a,b)仅仅表示一个方向，为了计算方便，我们取单位向量，即满足下式：
a2+b2=1a^2+b^2=1 a2+b2=1
第iii个点到直线的距离，可以表示为向量(xi−x0,yi−y0)(x_i-x_0,y_i-y_0)(xi−x0,yi−y0)在法向量方向(a,b)(a,b)(a,b)投影的长度，目标函数可以写成：
J2=∑indi2=∑in((xi−x0,yi−y0)⋅(a,b))2a2+b2=∑in[a(xi−x0)+b(yi−y0)]2J_2=\sum_{i}^{n}d_i^{2}=\sum_{i}^{n}\frac{((x_i-x_0,y_i-y_0)\cdot(a,b))^2}{a^2+b^2}=\sum_{i}^{n}[a(x_i-x_0)+b(y_i-y_0)]^2 J2=i∑ndi2=i∑na2+b2((xi−x0,yi−y0)⋅(a,b))2=i∑n[a(xi−x0)+b(yi−y0)]2
由极值点和可导函数之间的关系得，极值点可以推出一阶导数为0，那么一阶导数为0就是极值点的必要条件。那么可以推出：一阶导数不为0，不是极值点。

将目标函数J2J_2J2对x0x_0x0和y0y_0y0求偏导，并令其等于0，得：
∂J2∂x0=−2a∑in[a(xi−x0)+b(yi−y0)]=0\frac{\partial J_2}{\partial x_0}=-2a\sum_{i}^{n}[a(x_i-x_0)+b(y_i-y_0)]=0 ∂x0∂J2=−2ai∑n[a(xi−x0)+b(yi−y0)]=0

∂J2∂y0=−2b∑in[a(xi−x0)+b(yi−y0)]=0\frac{\partial J_2}{\partial y_0}=-2b\sum_{i}^{n}[a(x_i-x_0)+b(y_i-y_0)]=0 ∂y0∂J2=−2bi∑n[a(xi−x0)+b(yi−y0)]=0

上式等号两边同时除以n得
a(xˉ−x0)+b(yˉ−y0)=0a(\bar{x}-x_0)+b(\bar{y}-y_0)=0 a(xˉ−x0)+b(yˉ−y0)=0
可以看到点(xˉ,yˉ)(\bar{x},\bar{y})(xˉ,yˉ)满足直线方程，所以令x0=xˉ,y0=yˉx_0=\bar{x},y_0=\bar{y}x0=xˉ,y0=yˉ,此时目标函数变为：
J2=∑in[a(xi−xˉ)+b(yi−yˉ)]2=[ab][∑in(xi−xˉ)2∑in(xi−xˉ)(yi−yˉ)∑in(xi−xˉ)(yi−yˉ)∑in(yi−yˉ)2][ab]J_2=\sum_{i}^{n}[a(x_i-\bar{x})+b(y_i-\bar{y})]^2=\\ \begin{bmatrix} a & b \end{bmatrix} \begin{bmatrix} \sum_{i}^{n}(x_i-\bar{x})^2 &\sum_{i}^{n}(x_i-\bar{x})(y_i-\bar{y}) \\ \sum_{i}^{n}(x_i-\bar{x})(y_i-\bar{y}) &\sum_{i}^{n}(y_i-\bar{y})2 \end{bmatrix} \begin{bmatrix} a\\ b \end{bmatrix} J2=i∑n[a(xi−xˉ)+b(yi−yˉ)]2=[ab][∑in(xi−xˉ)2∑in(xi−xˉ)(yi−yˉ)∑in(xi−xˉ)(yi−yˉ)∑in(yi−yˉ)2][ab]
对目标函数J2J_2J2除以n可得：
J2n=[ab][SxxSxySxySyy][ab]=vTSv\frac{J_2}{n}= \begin{bmatrix} a & b \end{bmatrix} \begin{bmatrix} S_{xx} & S_{xy}\\ S_{xy} & S_{yy} \end{bmatrix} \begin{bmatrix} a\\ b \end{bmatrix}=v^TSv nJ2=[ab][SxxSxySxySyy][ab]=vTSv
其中SxxS_{xx}Sxx和SyyS_{yy}Syy为xxx和yyy的方差,SxyS_{xy}Sxy为xxx和yyy的协方差。

很明显，这是一个二次型求最小值的问题。相关定理证明也可参考线性代数及其应用 (豆瓣) (douban.com)

因为SSS是实对称矩阵，所以可以将其进行正交对角化分解：
S=[q1q2][λ100λ2][q1Tq2T]=QΛQTS= \begin{bmatrix} q_1 & q_2 \end{bmatrix} \begin{bmatrix} \lambda_1 &0 \\ 0 &\lambda_2 \end{bmatrix} \begin{bmatrix} q_{1}^{T}\\ q_{2}^{T} \end{bmatrix}=Q\Lambda Q^T S=[q1q2][λ100λ2][q1Tq2T]=QΛQT
有
J2n=vTSv=(vTQ)Λ(vTQ)T\frac{J_2}{n}=v^TSv=(v^TQ)\Lambda(v^TQ)^T nJ2=vTSv=(vTQ)Λ(vTQ)T
令u1=vTq1,u2=vTq2,u=[u1u2]u_1=v^Tq_1,u_2=v^Tq_2,u=\begin{bmatrix} u_1\\ u_2 \end{bmatrix}u1=vTq1,u2=vTq2,u=[u1u2]则：
J2n=uTΛu=λ1u12+λ2u22\frac{J_2}{n}=u^T\Lambda u=\lambda_1u_1^2+\lambda_2u_2^2 nJ2=uTΛu=λ1u12+λ2u22
因为：
uTu=vTQQTv=vTv=1u^Tu=v^TQQ^Tv=v^Tv=1 uTu=vTQQTv=vTv=1
所以uuu为单位矩阵，即u12+u22=1u_1^2+u_2^2=1u12+u22=1

不妨设λ1≤λ2\lambda_1\le\lambda_2λ1≤λ2,即当u1=1,u2=0u_1=1,u_2=0u1=1,u2=0时，J2J_2J2取得最小值λ1\lambda_1λ1，即v=q1v=q_1v=q1

所以最终的直线拟合结果是法向量vvv对应矩阵SSS的最小特征值的特征向量。

结果整理：
x0=xˉ,y0=yˉx_0=\bar{x},y_0=\bar{y} x0=xˉ,y0=yˉ
拟合直线法向量：
v=[ab]=q1v=\begin{bmatrix} a\\ b \end{bmatrix}=q_1 v=[ab]=q1
那么如何求矩阵SSS的最小特征值对应的特征向量呢

下面介绍使用奇异值分解求解矩阵SSS的最小特征值对应的特征向量，参考奇异值分解(SVD)原理与在降维中的应用 - 刘建平Pinard - 博客园 (cnblogs.com)

假设现在有nnn个二维点(x1,y1),(x2,y2)...(xn,yn)(x_1,y_1),(x_2,y_2)...(x_n,y_n)(x1,y1),(x2,y2)...(xn,yn),那么我们可以直接得到
A=[x1−xˉy1−yˉx2−xˉy2−yˉ......xn−xˉyn−yˉ]A=\begin{bmatrix} x_1-\bar{x}&y_1-\bar{y} \\ x_2-\bar{x}&y_2-\bar{y} \\ ...&... \\ x_n-\bar{x}&y_n-\bar{y} \end{bmatrix} A=⎣⎢⎢⎡x1−xˉx2−xˉ...xn−xˉy1−yˉy2−yˉ...yn−yˉ⎦⎥⎥⎤
那么我们可以观察得到，矩阵ATAA^TAATA就是矩阵SSS，那么根据奇异值分解的定义
A=UΣVTA=U\Sigma V^T A=UΣVT
上式中的VVV就是矩阵ATAA^TAATA也就是SSS的特征值对应的特征向量构成的正交矩阵，又因为在奇异值分解时，矩阵VVV中的特征值已经根据大小从大到小重新排列了，所以我们需要的就是 VVV的最后一列 ，也就是 VTV^TVT的最后一行 。

在python中，我们可以很方便的使用numpy中的np.linalg.svd求得相应的特征矩阵和特征向量。

从解齐次方程的角度推导

参考SVD之最小二乘【推导与证明】 - Brad_Lucas - 博客园 (cnblogs.com)和奇异值分解（SVD）方法求解最小二乘问题的原理 - 一抹烟霞 - 博客园 (cnblogs.com)

求解直线拟合的问题，可以转化为求解齐次方程Ax=0Ax=0Ax=0的问题，其中
A=[x1−xˉy1−yˉx2−xˉy2−yˉ......xn−xˉyn−yˉ]x=[ab]A=\begin{bmatrix} x_1-\bar{x}&y_1-\bar{y} \\ x_2-\bar{x}&y_2-\bar{y} \\ ...&... \\ x_n-\bar{x}&y_n-\bar{y} \end{bmatrix}\\x= \begin{bmatrix} a\\b \end{bmatrix} A=⎣⎢⎢⎡x1−xˉx2−xˉ...xn−xˉy1−yˉy2−yˉ...yn−yˉ⎦⎥⎥⎤x=[ab]
可以想到，当直线拟合的越好时，∣∣Ax∣∣||Ax||∣∣Ax∣∣越小

方法一

求min∣∣Ax∣∣||Ax||∣∣Ax∣∣，根据正交矩阵的性质：
∣∣Ax∣∣=∣∣UΣVTx∣∣=∣∣ΣVTx∣∣||Ax||=||U\Sigma V^Tx||=||\Sigma V^Tx|| ∣∣Ax∣∣=∣∣UΣVTx∣∣=∣∣ΣVTx∣∣
令y=VTxy=V^Txy=VTx，则：
∣∣ΣVTx∣∣=∣∣Σy∣∣||\Sigma V^Tx||=||\Sigma y|| ∣∣ΣVTx∣∣=∣∣Σy∣∣
进一步处理，求min(yTDTDy)min(y^TD^TDy)min(yTDTDy)
yTΣTΣy=σ12y12+⋯+σn2yn2σ1≥σ2≥…σn≥0y^T\Sigma^T\Sigma y=\sigma_1^2y_1^2+\dots+\sigma_n^2y_n^2\\ \sigma_1\ge\sigma_2\ge\dots\sigma_n\ge0 yTΣTΣy=σ12y12+⋯+σn2yn2σ1≥σ2≥…σn≥0
约束条件为∣∣y∣∣=1||y||=1∣∣y∣∣=1,如上所示，y=[00…1]y=\begin{bmatrix} 0&0&\dots&1\end{bmatrix}y=[00…1]是最小解，即：
x=Vy=[V1V2…Vn][000…1]=Vnx=Vy=\begin{bmatrix}V_1&V_2&\dots&V_n\end{bmatrix}\begin{bmatrix}0\\0\\0\\ \dots\\1\end{bmatrix}=V_n x=Vy=[V1V2…Vn]⎣⎢⎢⎢⎢⎡000…1⎦⎥⎥⎥⎥⎤=Vn
即最优解是VVV的最后一列

方法二

求min∣∣Ax∣∣=min(xTATAx)min||Ax||=min(x^TA^TAx)min∣∣Ax∣∣=min(xTATAx),其中A=UΣVT,UTU=I,VTV=IA=U\Sigma V^T,U^TU=I,V^TV=IA=UΣVT,UTU=I,VTV=I,,,Σ\SigmaΣ为对角阵。
ATA=VΣTUTUΣVT=VΣTΣVT=V[σmax2⋱σmin2]VTA^TA=V\Sigma^TU^TU\Sigma V^T=V\Sigma^T\Sigma V^T\\=V\begin{bmatrix}\sigma_{max}^2&&&\\ &\ddots&&\\&&&\sigma_{min}^2 \end{bmatrix}V^T\\ ATA=VΣTUTUΣVT=VΣTΣVT=V⎣⎡σmax2⋱σmin2⎦⎤VT
令V=[V0V1…Vn]V=\begin{bmatrix}V_0&V_1&\dots&V_n\end{bmatrix}V=[V0V1…Vn],得：
ATA=σmax2V0V0T+⋯+σmin2VnVnTA^TA=\sigma_{max}^2V_0V_0^T+\dots+\sigma_{min}^2V_nV_n^T ATA=σmax2V0V0T+⋯+σmin2VnVnT
由正交矩阵的性质：
VTV=IViTVj=0.i≠j∣∣Vi∣∣=1V^TV=I\\ V_i^TV_j=0.i\ne j\\ ||V_i||=1 VTV=IViTVj=0.i=j∣∣Vi∣∣=1
令x=∑0nkiVix=\sum_0^nk_iV_ix=∑0nkiVi，其中ViV_iVi为基向量
xTATAx=σmax2k02V0TV0V0TV0+⋯+σmin2kn2VnTVnVnTVn=k02σmax2+⋯+kn2σmin2x^TA^TAx=\sigma_{max}^2k_0^2V_0^TV_0V_0^TV_0+\dots+\sigma_{min}^2k_n^2V_n^TV_nV_n^TV_n\\=k_0^2\sigma_{max}^2+\dots+k_n^2\sigma_{min}^2 xTATAx=σmax2k02V0TV0V0TV0+⋯+σmin2kn2VnTVnVnTVn=k02σmax2+⋯+kn2σmin2
由于∣∣x∣∣=1||x||=1∣∣x∣∣=1，则
∑0nki2=1\sum_0^nk_i^2=1 0∑nki2=1
上式取最小的情况为kn=1,ki=0(i≠n)k_n=1,k_i=0(i\ne n)kn=1,ki=0(i=n)

此时x=Vnx=V_nx=Vn，为σmin\sigma_{min}σmin对应的列向量。

总结

我们发现，如果把上述的推导过程变为三维，把二维直线的点法式方程变为空间平面的点法式方程，上述推导依旧是成立的，所以我们已经得到了拟合空间平面的理论基础了。

References

RANSAC算法详解(附Python拟合直线模型代码) - 知乎 (zhihu.com)

线性拟合1-最小二乘 - 知乎 (zhihu.com)

线性拟合2-正交回归 - 知乎 (zhihu.com)

线性代数及其应用 (豆瓣) (douban.com)

奇异值分解(SVD)原理与在降维中的应用 - 刘建平Pinard - 博客园 (cnblogs.com)

奇异值分解（SVD）方法求解最小二乘问题的原理 - 一抹烟霞 - 博客园 (cnblogs.com)

SVD之最小二乘【推导与证明】 - Brad_Lucas - 博客园 (cnblogs.com)