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 + 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 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)进行整理得到。

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

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

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

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

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

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

∂ J 2 ∂ y 0 = − 2 b ∑ i n [ a ( x i − x 0 ) + b ( y i − y 0 ) ] = 0 frac{partial J_2}{partial y_0}=-2bsum_{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 ˉ − x 0 ) + b ( y ˉ − y 0 ) = 0 a(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ˉ​)满足直线方程，所以令 x 0 = x ˉ , y 0 = y ˉ x_0=bar{x},y_0=bar{y} x0​=xˉ,y0​=yˉ​,此时目标函数变为：
J 2 = ∑ i n [ a ( x i − x ˉ ) + b ( y i − y ˉ ) ] 2 = [ a b ] [ ∑ i n ( x i − x ˉ ) 2 ∑ i n ( x i − x ˉ ) ( y i − y ˉ ) ∑ i n ( x i − x ˉ ) ( y i − y ˉ ) ∑ i n ( y i − y ˉ ) 2 ] [ a b ] 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=[a​b​][∑in​(xi​−xˉ)2∑in​(xi​−xˉ)(yi​−yˉ​)​∑in​(xi​−xˉ)(yi​−yˉ​)∑in​(yi​−yˉ​)2​][ab​]
对目标函数 J 2 J_2 J2​除以n可得：
J 2 n = [ a b ] [ S x x S x y S x y S y y ] [ a b ] = v T S v 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​​=[a​b​][Sxx​Sxy​​Sxy​Syy​​][ab​]=vTSv
其中 S x x S_{xx} Sxx​和 S y y S_{yy} Syy​为 x x x和 y y y的方差, S x y S_{xy} Sxy​为 x x x和 y y y的协方差。

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

因为 S S S是实对称矩阵，所以可以将其进行正交对角化分解：
S = [ q 1 q 2 ] [ λ 1 0 0 λ 2 ] [ q 1 T q 2 T ] = Q Λ Q T S= 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}=QLambda Q^T S=[q1​​q2​​][λ1​0​0λ2​​][q1T​q2T​​]=QΛQT
有
J 2 n = v T S v = ( v T Q ) Λ ( v T Q ) T frac{J_2}{n}=v^TSv=(v^TQ)Lambda(v^TQ)^T nJ2​​=vTSv=(vTQ)Λ(vTQ)T
令 u 1 = v T q 1 , u 2 = v T q 2 , u = [ u 1 u 2 ] u_1=v^Tq_1,u_2=v^Tq_2,u=begin{bmatrix} u_1\ u_2 end{bmatrix} u1​=vTq1​,u2​=vTq2​,u=[u1​u2​​]则：
J 2 n = u T Λ u = λ 1 u 1 2 + λ 2 u 2 2 frac{J_2}{n}=u^TLambda u=lambda_1u_1^2+lambda_2u_2^2 nJ2​​=uTΛu=λ1​u12​+λ2​u22​
因为：
u T u = v T Q Q T v = v T v = 1 u^Tu=v^TQQ^Tv=v^Tv=1 uTu=vTQQTv=vTv=1
所以 u u u为单位矩阵，即 u 1 2 + u 2 2 = 1 u_1^2+u_2^2=1 u12​+u22​=1

不妨设 λ 1 ≤ λ 2 lambda_1lelambda_2 λ1​≤λ2​,即当 u 1 = 1 , u 2 = 0 u_1=1,u_2=0 u1​=1,u2​=0时， J 2 J_2 J2​取得最小值 λ 1 lambda_1 λ1​，即 v = q 1 v=q_1 v=q1​

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

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

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

假设现在有 n n n个二维点 ( x 1 , y 1 ) , ( x 2 , y 2 ) . . . ( x n , y n ) (x_1,y_1),(x_2,y_2)...(x_n,y_n) (x1​,y1​),(x2​,y2​)...(xn​,yn​),那么我们可以直接得到
A = [ x 1 − x ˉ y 1 − y ˉ x 2 − x ˉ y 2 − y ˉ . . . . . . x n − x ˉ y n − 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ˉ​​⎦⎥⎥⎤​
那么我们可以观察得到，矩阵 A T A A^TA ATA就是矩阵 S S S，那么根据奇异值分解的定义
A = U Σ V T A=USigma V^T A=UΣVT
上式中的 V V V就是矩阵 A T A A^TA ATA也就是 S S S的特征值对应的特征向量构成的正交矩阵，又因为在奇异值分解时，矩阵 V V V中的特征值已经根据大小从大到小重新排列了，所以我们需要的就是** V V V的最后一列**，也就是** V T V^T VT的最后一行**。

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

从解齐次方程的角度推导

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

求解直线拟合的问题，可以转化为求解齐次方程 A x = 0 Ax=0 Ax=0的问题，其中
A = [ x 1 − x ˉ y 1 − y ˉ x 2 − x ˉ y 2 − y ˉ . . . . . . x n − x ˉ y n − y ˉ ] x = [ a b ] 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​]
可以想到，当直线拟合的越好时，$left | Ax right | $越小

方法一

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

方法二

求 m i n ∣ ∣ A x ∣ ∣ = m i n ( x T A T A x ) min||Ax||=min(x^TA^TAx) min∣∣Ax∣∣=min(xTATAx),其中 A = U Σ V T , U T U = I , V T V = I A=USigma V^T,U^TU=I,V^TV=I A=UΣVT,UTU=I,VTV=I , , , Σ Sigma Σ为对角阵。
A T A = V Σ T U T U Σ V T = V Σ T Σ V T = V [ σ m a x 2 ⋱ σ m i n 2 ] V T A^TA=VSigma^TU^TUSigma V^T=VSigma^TSigma V^T\=Vbegin{bmatrix}sigma_{max}^2&&&\ &ddots&&\&&&sigma_{min}^2 end{bmatrix}V^T\ ATA=VΣTUTUΣVT=VΣTΣVT=V⎣⎡​σmax2​​⋱​​σmin2​​⎦⎤​VT
令 V = [ V 0 V 1 … V n ] V=begin{bmatrix}V_0&V_1&dots&V_nend{bmatrix} V=[V0​​V1​​…​Vn​​],得：
A T A = σ m a x 2 V 0 V 0 T + ⋯ + σ m i n 2 V n V n T A^TA=sigma_{max}^2V_0V_0^T+dots+sigma_{min}^2V_nV_n^T ATA=σmax2​V0​V0T​+⋯+σmin2​Vn​VnT​
由正交矩阵的性质：
V T V = I V i T V j = 0. i ≠ j ∣ ∣ V i ∣ ∣ = 1 V^TV=I\ V_i^TV_j=0.ine j\ ||V_i||=1 VTV=IViT​Vj​=0.i​=j∣∣Vi​∣∣=1
令 x = ∑ 0 n k i V i x=sum_0^nk_iV_i x=∑0n​ki​Vi​，其中 V i V_i Vi​为基向量
x T A T A x = σ m a x 2 k 0 2 V 0 T V 0 V 0 T V 0 + ⋯ + σ m i n 2 k n 2 V n T V n V n T V n = k 0 2 σ m a x 2 + ⋯ + k n 2 σ m i n 2 x^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^2sigma_{max}^2+dots+k_n^2sigma_{min}^2 xTATAx=σmax2​k02​V0T​V0​V0T​V0​+⋯+σmin2​kn2​VnT​Vn​VnT​Vn​=k02​σmax2​+⋯+kn2​σmin2​
由于 ∣ ∣ x ∣ ∣ = 1 ||x||=1 ∣∣x∣∣=1，则
∑ 0 n k i 2 = 1 sum_0^nk_i^2=1 0∑n​ki2​=1
上式取最小的情况为 k n = 1 , k i = 0 ( i ≠ n ) k_n=1,k_i=0(ine n) kn​=1,ki​=0(i​=n)

此时 x = V n x=V_n x=Vn​，为 σ m i n 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)