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

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

$$y = k x + b$$

那么

$$k= \frac{\sum{(x_i-\bar x)(y_i-\bar y)}}{\sum{(x_i - \bar x)^2}}$$
$$b = \bar y - k \bar x$$

这个关系式许多教科书上都有详细的推导，无需多说。
今天要说的是另一种情况，当我们的数据有可能落在一条竖直的直线上，也就是$k$ 有可能为$\infty$ 时，应该如何做拟合。这时我们肯定就不能用$y = k x + b$ 了，但是可以将这个表达式变变形。我们知道

$$k = \tan\theta = \frac{\sin \theta}{\cos \theta}$$
那么原来的直线方程可以写为：
$$y \cos \theta = x \sin \theta + b \cos \theta$$
或者写为更一般的形式：
$$a x + b y + c = 0$$
同时满足附加条件：
$$a^2 + b^2 = 1$$
下面就来说说这种形式的直线方程如何拟合。

点到直线的垂直距离

首先先要解决一个小问题，一个点 $(x_i, y_i)$ 到这条直线的距离是多少。
直接计算有点麻烦，我们先考虑一种简单的形式，点 $(0, 0)$ 到这条直线的距离。

直线方程：

$$x \sin(\theta) - y \cos(\theta) + c = 0$$

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

$$x \sin(\theta + \pi / 2) - y \cos(\theta + \pi / 2) = 0$$

化简后得到：

$$x \cos(\theta) + y \sin(\theta) = 0$$

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

$$\begin{cases} x \sin(\theta) - y \cos(\theta) + c = 0 \x \cos(\theta) + y \sin(\theta) = 0 \end{cases}$$

简单计算后可以得到：

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

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

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

$$\begin{cases} x‘ = x - x_i \y‘ = y - y_i \end{cases}$$

那么原始坐标系下的点$(x_i, y_i)$ 在这个新的坐标系下称为了坐标原点$(0,0)$。
原始坐标系下的直线方程 $a x + b y + c = 0$ 成为了：

$$a (x‘ + x_i) + b (y‘ + y_i) + c = 0$$

整理一下成为：

$$a x‘ + b y‘ + (a x_i + b y_i + c) = 0$$

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

直线拟合

我们要求的直线方程的系数 $a$、$b$ 和 $c$ 就是在

$$a^2 + b^2 = 1$$ 的条件下，使得下面求和式：
$$f = \sum (a x_i + b y_i + c)^2$$
取得极小值的 $a$、$b$ 和 $c$。

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

$$f = \sum (a x_i + b y_i + c)^2 - \lambda (a^2 + b^2 - 1)$$

$$\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}$$

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

$$\frac{\partial f}{\partial c} = 2 \times \sum (a x_i + b y_i + c) = 0$$

变形之后可以得到：

$$a \bar x + b \bar y + c =0 \\bar x = \sum x_i / N \\bar y = \sum y_i / N$$
用这个式子， $c$ 就确定了。可以将确定后的 $c$ 带入 $f$ 的表达式。

$$f = \sum (a (x_i - \bar x) + b (y_i - \bar y))^2 - \lambda (x^2 + y^2 - 1)$$

$$\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}$$

$$\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}$$

其中：

$$\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}$$

所以，我们的 $a$、$b$ 和 $\lambda$ 应该满足如下线性方程组。

$$\begin{cases} a D_{xx} + b D_{xy} = a \lambda \a D_{xy} + b D_{yy} = b \lambda \end{cases}$$

或者写成矩阵形式：

$$\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$ 是 $\left ( \begin{array}{ccc} D_{xx} & D_{xy} \D_{xy} & D_{yy} \end{array} \right )$ 的特征值， $a , b$ 是 $\left ( \begin{array}{ccc} D_{xx} & D_{xy} \D_{xy} & D_{yy} \end{array} \right )$
的特征向量。

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

$$\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$$

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

$$\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$$
有非零解的条件是：
$$\left | \begin{array}{ccc} D_{xx} -\lambda & D_{xy} \D_{xy} & D_{yy} -\lambda \end{array} \right | = 0$$

$$\lambda = \frac{(D_{xx} + D_{yy}) \pm \sqrt{(D_{xx} - D_{yy})^2 + 4 D_{xy}^2} }{2}$$

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

$$\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 )$$

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

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

后记

写到这里关于直线拟合算法的实现方法就算是说清楚了。但其实还有一个重要的问题没有交代。我们做一元线性回归时会去求相关系数。相关系数反映的是数据点与拟合直线的吻合程度，对于我们的算法，也应该提出一个类似的指标。另外，求出的系数 a、b、c 也应给出不确定度或者置信区间。
关于这些问题，我准备再开一篇博客来详细的说说。（其实是因为有些问题还没想好如何解决(^.^)）