RANSAC算法与原理（二）

之前介绍了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 + 10

fig = 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，b
        a = (y_2 - y_1) / (x_2 - x_1)
        b = y_1 - a * x_1

        # 算出内点数目
        total_inlier = 0
        for index in range(SIZE * 2):
            y_estimate = a * RANDOM_X[index] + b
            if 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_inlier
            best_a = a
            best_b = b

        # 判断是否当前模型已经符合超过一半的点
        if total_inlier > SIZE:
            break
	else:
        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] + b
            if 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)进行整理得到。

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

横坐标和纵坐标的误差定义如下：
$$\widehat{{\eta }_i}=x_i-\widehat{x_i}$$

$$\widehat{{ϵ}_i}=y_i-\widehat{y_i}$$

综合考虑横纵坐标的误差，目标函数的形式如下：
$$J_2=\sum _i^n[(\widehat{{\eta }_i}{)}^2+(\widehat{{ϵ}_i}{)}^2]=\sum _i^n[(x_i-\widehat{x_i}{)}^2+(y_i-\widehat{y_i}{)}^2]$$
因为要求目标函数的最小值，所以就是求观测点$(x_i,y_i)$到估计点$(\widehat{x_i},\widehat{y_i})$的最小值，因为点$(x_i,y_i)$我们是已知的，也就是求点$(x_i,y_i)$到拟合直线的距离的和的最小值。估计点$(\widehat{x_i},\widehat{y_i})$是点$(x_i,y_i)$在拟合直线上的正交投影点。所以目标函数可以写成：
$$J_2=\sum _i^n{d}_i^2$$
$d_i$为第$i$个点$(x_i,y_i)$到拟合直线的距离。

用点法式拟合直线的表示形式为
$$a(x-x_0)+b(y-y_0)=0$$
其中，$(x_0,y_0)$为经过直线的一个点，$(a,b)$表示直线的法向量。$(a,b)$仅仅表示一个方向，为了计算方便，我们取单位向量，即满足下式：
$$a^2+b^2=1$$
第$i$个点到直线的距离，可以表示为向量$(x_i-x_0,y_i-y_0)$在法向量方向$(a,b)$投影的长度，目标函数可以写成：
$$J_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$$
由极值点和可导函数之间的关系得，极值点可以推出一阶导数为0，那么一阶导数为0就是极值点的必要条件。那么可以推出：一阶导数不为0，不是极值点。

将目标函数$J_2$对$x_0$和$y_0$求偏导，并令其等于0，得：
$$\frac{\partial J_2}{\partial x_0}=-2a\sum _i^n[a(x_i-x_0)+b(y_i-y_0)]=0$$

$$\frac{\partial J_2}{\partial y_0}=-2b\sum _i^n[a(x_i-x_0)+b(y_i-y_0)]=0$$

上式等号两边同时除以n得
$$a(\bar{x}-x_0)+b(\bar{y}-y_0)=0$$
可以看到点$(\bar{x},\bar{y})$满足直线方程，所以令$x_0=\bar{x},y_0=\bar{y}$,此时目标函数变为：
$$\begin{gathered} J_2=\sum _i^n[a(x_i-\bar{x})+b(y_i-\bar{y}){]}^2= \\ \left[\begin{array}{cc}a& b\end{array}\right]\left[\begin{array}{cc}{\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{array}\right]\left[\begin{array}{c}a\\ b\end{array}\right] \end{gathered}$$
对目标函数$J_2$除以n可得：
$$\frac{J_2}{n}=\left[\begin{array}{cc}a& b\end{array}\right]\left[\begin{array}{cc}{S}_{xx}& S_{xy}\\ S_{xy}& S_{yy}\end{array}\right]\left[\begin{array}{c}a\\ b\end{array}\right]=v^{T}Sv$$
其中$S_{xx}$和$S_{yy}$为$x$和$y$的方差,$S_{xy}$为$x$和$y$的协方差。

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

因为$S$是实对称矩阵，所以可以将其进行正交对角化分解：
$$S=\left[\begin{array}{cc}{q}_1& q_2\end{array}\right]\left[\begin{array}{cc}{\lambda }_1& 0\\ 0& {\lambda }_2\end{array}\right]\left[\begin{array}{c}{q}_1^T\\ q_2^T\end{array}\right]=Q\Lambda Q^T$$
有
$$\frac{J_2}{n}=v^{T}Sv=(v^{T}Q)\Lambda (v^{T}Q{)}^T$$
令$u_1=v^T{q}_1,u_2=v^T{q}_2,u=\left[\begin{array}{c}{u}_1\\ u_2\end{array}\right]$则：
$$\frac{J_2}{n}=u^T\Lambda u={\lambda }_1{u}_1^2+{\lambda }_2{u}_2^2$$
因为：
$$u^{T}u=v^{T}Q{Q}^{T}v=v^{T}v=1$$
所以$u$为单位矩阵，即$u_1^2+u_2^2=1$

不妨设${\lambda }_1\le {\lambda }_2$,即当$u_1=1,u_2=0$时，$J_2$取得最小值${\lambda }_1$，即$v=q_1$

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

结果整理：
$$x_0=\bar{x},y_0=\bar{y}$$
拟合直线法向量：
$$v=\left[\begin{array}{c}a\\ b\end{array}\right]=q_1$$
那么如何求矩阵$S$的最小特征值对应的特征向量呢

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

假设现在有$n$个二维点$(x_1,y_1),(x_2,y_2)...(x_n,y_n)$,那么我们可以直接得到
$$A=\left[\begin{array}{cc}{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{array}\right]$$
那么我们可以观察得到，矩阵$A^{T}A$就是矩阵$S$，那么根据奇异值分解的定义
$$A=U\Sigma V^T$$
上式中的$V$就是矩阵$A^{T}A$也就是$S$的特征值对应的特征向量构成的正交矩阵，又因为在奇异值分解时，矩阵$V$中的特征值已经根据大小从大到小重新排列了，所以我们需要的就是 $V$的最后一列 ，也就是 $V^T$的最后一行 。

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

从解齐次方程的角度推导

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

求解直线拟合的问题，可以转化为求解齐次方程$Ax=0$的问题，其中
$$\begin{gathered} A=\left[\begin{array}{cc}{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{array}\right] \\ x=\left[\begin{array}{c}a\\ b\end{array}\right] \end{gathered}$$
可以想到，当直线拟合的越好时，$||Ax||$越小

方法一

求min$||Ax||$，根据正交矩阵的性质：
$$||Ax||=||U\Sigma V^{T}x||=||\Sigma V^{T}x||$$
令$y=V^{T}x$，则：
$$||\Sigma V^{T}x||=||\Sigma y||$$
进一步处理，求$min(y^T{D}^{T}Dy)$
$$\begin{gathered} y^T{\Sigma }^T\Sigma y={\sigma }_1^2{y}_1^2+\cdots +{\sigma }_n^2{y}_n^2 \\ {\sigma }_1\ge {\sigma }_2\ge \dots {\sigma }_n\ge 0 \end{gathered}$$
约束条件为$||y||=1$,如上所示，$y=\left[\begin{array}{cccc}0& 0& \dots & 1\end{array}\right]$是最小解，即：
$$x=Vy=\left[\begin{array}{cccc}{V}_1& V_2& \dots & V_n\end{array}\right]\left[\begin{array}{c}0\\ 0\\ 0\\ \dots \\ 1\end{array}\right]=V_n$$
即最优解是$V$的最后一列

方法二

求$min||Ax||=min(x^T{A}^{T}Ax)$,其中$A=U\Sigma V^T,U^{T}U=I,V^{T}V=I$$,$$\Sigma$为对角阵。
$$\begin{gathered} A^{T}A=V{\Sigma }^T{U}^{T}U\Sigma V^T=V{\Sigma }^T\Sigma V^T \\ =V\left[\begin{array}{cccc}{\sigma }_{max}^2& & & \\ & \ddots & & \\ & & & {\sigma }_{min}^2\end{array}\right]V^T \end{gathered}$$
令$V=\left[\begin{array}{cccc}{V}_0& V_1& \dots & V_n\end{array}\right]$,得：
$$A^{T}A={\sigma }_{max}^2{V}_0{V}_0^T+\cdots +{\sigma }_{min}^2{V}_n{V}_n^T$$
由正交矩阵的性质：
$$\begin{gathered} V^{T}V=I \\ {V}_i^T{V}_j=0.i\ne j \\ ||V_i||=1 \end{gathered}$$
令$x={\sum }_0^n{k}_i{V}_i$，其中$V_i$为基向量
$$\begin{gathered} x^T{A}^{T}Ax={\sigma }_{max}^2{k}_0^2{V}_0^T{V}_0{V}_0^T{V}_0+\cdots +{\sigma }_{min}^2{k}_n^2{V}_n^T{V}_n{V}_n^T{V}_n \\ =k_0^2{\sigma }_{max}^2+\cdots +k_n^2{\sigma }_{min}^2 \end{gathered}$$
由于$||x||=1$，则
$$\sum _0^n{k}_i^2=1$$
上式取最小的情况为$k_n=1,k_i=0(i\ne n)$

此时$x=V_n$，为${\sigma }_{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)