ax+by+c=0 型直线拟合算法

所谓直线拟合，通常也叫做线性拟合、一元线性回归。指的是当我们有一批数据(xi,yi)(x_i, y_i)，这些数据在平面坐标系下落在一条直线上，或近似的落在一条直线上。我们就要求出这条直线的参数。如果这条直线可以写为：y=kx+by = k x + b

那么

k=∑(xi−x¯)(yi−y¯)∑(xi−x¯)2k= \frac{\sum{(x_i-\bar x)(y_i-\bar y)}}{\sum{(x_i - \bar x)^2}}


b=y¯−kx¯b = \bar y - k \bar x

这个关系式许多教科书上都有详细的推导，无需多说。

今天要说的是另一种情况，当我们的数据有可能落在一条竖直的直线上，也就是kk 有可能为∞\infty 时，应该如何做拟合。这时我们肯定就不能用y=kx+by = k x + b 了，但是可以将这个表达式变变形。我们知道

k=tanθ=sinθcosθk = \tan\theta = \frac{\sin \theta}{\cos \theta}

那么原来的直线方程可以写为：

ycosθ=xsinθ+bcosθy \cos \theta = x \sin \theta + b \cos \theta

或者写为更一般的形式：

ax+by+c=0 a x + b y + c = 0

同时满足附加条件：

a2+b2=1 a^2 + b^2 = 1

下面就来说说这种形式的直线方程如何拟合。

点到直线的垂直距离

首先先要解决一个小问题，一个点 (xi,yi)(x_i, y_i) 到这条直线的距离是多少。

直接计算有点麻烦，我们先考虑一种简单的形式，点 (0,0)(0, 0) 到这条直线的距离。

直线方程：

xsin(θ)−ycos(θ)+c=0 x \sin(\theta) - y \cos(\theta) + c = 0

那么与它垂直的过原点的直线方程为：

xsin(θ+π/2)−ycos(θ+π/2)=0
x \sin(\theta + \pi / 2) - y \cos(\theta + \pi / 2) = 0

化简后得到：

xcos(θ)+ysin(θ)=0
x \cos(\theta) + y \sin(\theta) = 0

这两个直线的交点为两条直线方程组成的方程组的解：

{xsin(θ)−ycos(θ)+c=0xcos(θ)+ysin(θ)=0
\begin{cases}
x \sin(\theta) - y \cos(\theta) + c = 0 \\
x \cos(\theta) + y \sin(\theta) = 0
\end{cases}

简单计算后可以得到：

{x=−sin(θ)cy=cos(θ)c
\begin{cases}
x = - \sin(\theta) c \\
y = \cos(\theta) c
\end{cases}

容易看出交点到原点的距离为 cc，而这就是原点到直线的距离。

下面就可以求点 (xi,yi)(x_i, y_i) 到这条直线的距离了。只需要做个坐标变换。新坐标系下的坐标变量为 x′x' 和 y′y'。 新旧坐标系的关系如下。

{x′=x−xiy′=y−yi
\begin{cases}
x' = x - x_i \\
y' = y - y_i
\end{cases}

那么原始坐标系下的点(xi,yi)(x_i, y_i) 在这个新的坐标系下称为了坐标原点(0,0)(0,0)。

原始坐标系下的直线方程 ax+by+c=0a x + b y + c = 0 成为了：

a(x′+xi)+b(y′+yi)+c=0 a (x' + x_i) + b (y' + y_i) + c = 0

整理一下成为：

ax′+by′+(axi+byi+c)=0 a x' + b y' + (a x_i + b y_i + c) = 0

在这个坐标系下直线到原点的距离为： (axi+byi+c)(a x_i + b y_i + c)，这也就是我们要求的直线到点的垂直距离。

直线拟合

我们要求的直线方程的系数 aa、bb 和 cc 就是在

a2+b2=1 a^2 + b^2 = 1 的条件下，使得下面求和式：

f=∑(axi+byi+c)2 f = \sum (a x_i + b y_i + c)^2

取得极小值的 aa、bb 和 cc。

这个问题有标准的处理方法：拉格朗日乘子法。

f=∑(axi+byi+c)2−λ(a2+b2−1) f = \sum (a x_i + b y_i + c)^2 - \lambda (a^2 + b^2 - 1)

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪∂f∂a=0∂f∂b=0∂f∂c=0∂f∂λ=0
\begin{cases}
\frac{\partial f}{\partial a} = 0 \\
\frac{\partial f}{\partial b} = 0 \\
\frac{\partial f}{\partial c} = 0 \\
\frac{\partial f}{\partial {\lambda}} = 0
\end{cases}

这里先不急着把这几个式子都展开了。只来看看 ∂f∂c\frac{\partial f}{\partial c} 这一项。

∂f∂c=2×∑(axi+byi+c)=0
\frac{\partial f}{\partial c} = 2 \times \sum (a x_i + b y_i + c) = 0

变形之后可以得到：

ax¯+by¯+c=0x¯=∑xi/Ny¯=∑yi/N
a \bar x + b \bar y + c =0 \\
\bar x = \sum x_i / N \\
\bar y = \sum y_i / N


用这个式子， cc 就确定了。可以将确定后的 cc 带入 ff 的表达式。

f=∑(a(xi−x¯)+b(yi−y¯))2−λ(x2+y2−1)
f = \sum (a (x_i - \bar x) + b (y_i - \bar y))^2 - \lambda (x^2 + y^2 - 1)

∂f∂a=2∑(a(xi−x¯)+b(yi−y¯))(xi−x¯)−2aλ=2(aDxx+bDxy)−2aλ=0
\begin{split}
\frac{\partial f}{\partial a} &= 2 \sum (a (x_i - \bar x) + b (y_i - \bar y)) (x_i - \bar x) - 2 a \lambda \\
&= 2 (a D_{xx} + b D_{xy}) - 2 a \lambda \\
& = 0
\end{split}

∂f∂b=2∑(a(xi−x¯)+b(yi−y¯))(yi−x¯)−2bλ=2(aDxy+bDyy)−2bλ=0
\begin{split}
\frac{\partial f}{\partial b} &= 2 \sum (a (x_i - \bar x) + b (y_i - \bar y)) (y_i - \bar x) - 2 b \lambda \\
&= 2 (a D_{xy} + b D_{yy}) - 2 b \lambda \\
& = 0
\end{split}

其中：

⎧⎩⎨⎪⎪Dxx=∑(xi−x¯)2Dxy=∑(xi−x¯)(yi−y¯)Dyy=∑(yi−y¯)2
\begin{cases}
D_{xx} = \sum (x_i - \bar x)^2 \\
D_{xy} = \sum (x_i - \bar x) ( y_i - \bar y) \\
D_{yy} = \sum (y_i - \bar y)^2
\end{cases}

所以，我们的 aa、bb 和 λ\lambda 应该满足如下线性方程组。

{aDxx+bDxy=aλaDxy+bDyy=bλ
\begin{cases}
a D_{xx} + b D_{xy} = a \lambda \\
a D_{xy} + b D_{yy} = b \lambda
\end{cases}

或者写成矩阵形式：

(DxxDxyDxyDyy)(ab)=λ(ab) \left (
\begin{array}{ccc}
D_{xx} & D_{xy} \\
D_{xy} & D_{yy}
\end{array} \right ) \left (
\begin{array}{ccc}
a \\
b
\end{array} \right )
= \lambda \left (
\begin{array}{ccc}
a \\
b
\end{array} \right )

这样写就很清楚了。λ\lambda 是 (DxxDxyDxyDyy)\left (
\begin{array}{ccc}
D_{xx} & D_{xy} \\
D_{xy} & D_{yy}
\end{array} \right ) 的特征值， a,ba , b 是 (DxxDxyDxyDyy)\left (
\begin{array}{ccc}
D_{xx} & D_{xy} \\
D_{xy} & D_{yy}
\end{array} \right )

的特征向量。

但是这个矩阵有两个特征值，我们应该选用哪个特征值呢。

∑(a(xi−x¯)+b(yi−y¯))2=(ab)(DxxDxyDxyDyy)(ab)=λ(a2+b2)=λ
\sum \left(a (x_i - \bar x) + b (y_i - \bar y) \right)^2
= \left (\begin{array}{ccc}
a & b
\end{array} \right )
\left ( \begin{array}{ccc}
D_{xx} & D_{xy} \\
D_{xy} & D_{yy}
\end{array} \right )
\left ( \begin{array}{ccc}
a \\
b
\end{array} \right ) = \lambda (a^2 + b^2) = \lambda

所以，我们应该选取两个特征值中较小的那个。

(Dxx−λDxyDxyDyy−λ)(ab)=0 \left (
\begin{array}{ccc}
D_{xx} -\lambda & D_{xy} \\
D_{xy} & D_{yy} -\lambda
\end{array} \right ) \left (
\begin{array}{ccc}
a \\
b
\end{array} \right )
= 0


有非零解的条件是：

∣∣∣Dxx−λDxyDxyDyy−λ∣∣∣=0 \left |
\begin{array}{ccc}
D_{xx} -\lambda & D_{xy} \\
D_{xy} & D_{yy} -\lambda
\end{array} \right |
= 0

λ=(Dxx+Dyy)±(Dxx−Dyy)2+4D2xy−−−−−−−−−−−−−−−−−√2
\lambda = \frac{(D_{xx} + D_{yy}) \pm \sqrt{(D_{xx} - D_{yy})^2 + 4 D_{xy}^2} }{2}

较小的那个特征根是： λ=(Dxx+Dyy)−(Dxx−Dyy)2+4D2xy√2\lambda = \frac{(D_{xx} + D_{yy}) - \sqrt{(D_{xx} - D_{yy})^2 + 4 D_{xy}^2} }{2}

(ab)=1D2xy+(λ−Dxx)2−−−−−−−−−−−−−−−√(Dxyλ−Dxx)
\left (\begin{array}{ccc}
a \\ b
\end{array} \right )
= \frac{1}{\sqrt{D_{xy}^2+(\lambda-D{xx})^2}}
\left (\begin{array}{ccc}
D_{xy} \\ \lambda - D_{xx}
\end{array} \right )

至此，这个直线拟合问题的解决了。但是，关于这个直线拟合与传统的一元线性回归算法的区别我还想说几句。

传统的一元线性回归是认为数据点 (xi,yi)(x_i, y_i) 中只有 yiy_i 是有误差的，因此确定点到直线距离时用的是 y 方向的距离。 我本文中的算法认为数据点 (xi,yi)(x_i, y_i) 都是有误差的，并且不确定度是相同的，因此，数据点到直线的距离用的是垂直距离。这两个异同可以参考下面的图片。当这条直线的斜率很小时，这两种方法求得直线方程很接近，当直线斜率很大时，两个结果可能有很大的区别。具体应该用哪种还是要根据数据点的性质来确定。



下面给一个 C++ 的测试代码，用到了 Qt 的一些功能。当然也可以很容易的去除 Qt 的依赖。

#include <QCoreApplication>
#include <QVector>
#include <QPoint>
#include <math.h>
#include <iostream>
using namespace std;
void lineFit(const QVector<QPoint> &points, double &a, double &b, double &c);

int main(int argc, char *argv[])
{
    QCoreApplication app(argc, argv);
    QVector<QPoint> p;
    for(int i = 0; i < 20; i++)
    {
        p.append(QPoint (i , 2 * i + 3));
    }

    double a, b, c;
    lineFit(p, a, b, c);
    c = c/ b;
    a = a/ b;
    b = b/ b;

    cout << "line: a = " << a << "b = " << b << "c = " << c;
    return app.exec();
}

void lineFit(const QVector<QPoint> &points, double &a, double &b, double &c)
{
    int size = points.size();
    if(size == 0)
    {
        a = 0;
        b = 0;
        c = 0;
        return;
    }
    double x_mean = 0;
    double y_mean = 0;
    for(int i = 0; i < size; i++)
    {
        x_mean += points[i].x();
        y_mean += points[i].y();
    }
    x_mean /= size;
    y_mean /= size; //至此，计算出了 x y 的均值

    double Dxx = 0, Dxy = 0, Dyy = 0;

    for(int i = 0; i < size; i++)
    {
        Dxx += (points[i].x() - x_mean) * (points[i].x() - x_mean);
        Dxy += (points[i].x() - x_mean) * (points[i].y() - y_mean);
        Dyy += (points[i].y() - y_mean) * (points[i].y() - y_mean);
    }
    double lambda = ( (Dxx + Dyy) - sqrt( (Dxx - Dyy) * (Dxx - Dyy) + 4 * Dxy * Dxy) ) / 2.0;
    double den = sqrt( Dxy * Dxy + (lambda - Dxx) * (lambda - Dxx) );
    a = Dxy / den;
    b = (lambda - Dxx) / den;
    c = - a * x_mean - b * y_mean;
}

后记

写到这里关于直线拟合算法的实现方法就算是说清楚了。但其实还有一个重要的问题没有交代。我们做一元线性回归时会去求相关系数。相关系数反映的是数据点与拟合直线的吻合程度，对于我们的算法，也应该提出一个类似的指标。另外，求出的系数 a、b、c 也应给出不确定度或者置信区间。

关于这些问题，我准备再开一篇博客来详细的说说。（其实是因为有些问题还没想好如何解决(^.^)）