最小二乘法直线拟合:Ax+By+C=0

因为\( y =kx+b\) 在斜率无穷大或接近无穷大时的数值计算问题,所以在直线方程的选择上选用更一般的形式:$$Ax+By+C=0$$
已知$$  (x_1,y_1),(x_2,y_2),…,(x_n,y_n)$$求未知数\(A,B,C\)使得 \(x_i\) 到 直线 $$Ax+By+C=0$$ 的距离的平方和最小
$$ f=\sum_{i=1}^n (\frac{|Ax_i+By_i+C|}{\sqrt{A^2+B^2}})^2$$因为A,B可以成比例的变化,所以可以添加一个约束条件:$$A^2+B^2=1$$所以上面算式可以简化为
$$f=\sum_{i=1}^n(Ax_i+By_i+C)^2$$令\(f\) 对 \(C\)求偏导数得到$$ \frac{\partial f}{\partial C}=2(A\sum x+B\sum y+nC)=0$$ 解\(C\)得到 \( C=-\frac{A\sum x +  B\sum y}{n} \)。即$$C=-A\bar x -B \bar y$$。带入上面算式
$$ f=\sum(A(x_i-\bar x)+B(y_i -\bar y))^2 = \sum (Ap+Bq)^2 $$ $$ (p_i = x_i-\bar x,q_i = y_i – \bar y)$$
可以写成 矩阵乘法形式
$$f = \begin{bmatrix}A\\B \end{bmatrix}^T\begin{bmatrix}\sum p^2& \sum pq\\ \sum pq & \sum q^2\end{bmatrix} \begin{bmatrix}A\\B\end{bmatrix}$$
对应简写为 \( f = v^TSv\)。写到这里,问题转变为二次型求最小值的问题。在线性代数课本上可能已经写出了问题的解,即 \( v=\begin{bmatrix} A\\B\end{bmatrix}\)的值为S对应最小特征值时的特征向量。
采用Eigen库的EigenSolver可以计算特征值和特征向量。值得提到的是 Eigen 的eigenvalues()</方法返回的是std::complex<double> 复数数组,因为 \( S\) 为实对称矩阵,所以其特征值为 实数,即取 complex<double>::real()方法就可以了。
解出了A,B ,参数C通过开始的公式可以求出。

在实际应用中,有小伙伴希望给出 为什么v的值是对应S特征值为最小时的特征向量的值,这里我尝试给出说明。
因为\(S\) 为实对称矩阵,其可以对角化分解为 \( S= PDP^{-1}=PDP^T\)$$ f=v^TSv = (v^TP)D(v^TP)^T$$
D为特征值对角矩阵,P为对应的特征向量。
\(v^tP =[ v^tP_1,v^tP_2] =[u_1,u_2]\)为向量v与p的点积,其范围在-1到1之间。
$$f= u_1^2\lambda _1 + u_2^2\lambda _2$$,因为\( p_0,p1\)是正交的,v的模长为1,可以假设\( p_1,p_2\)的坐标为\((1,0),(0,1)\),v的坐标为\((x,y),  x^2+y^2=1\) 。\(\lambda_1>\lambda_2\)

$$f=u_1^2\lambda _1 + u_2^2\lambda _2 = x^2\lambda _1 + y^2\lambda _2 = \lambda_2 + (\lambda_1-\lambda_2)x^2$$

当 \((x,y) = (0,1) \) 即当v与较小的\(\lambda\)对应的\(p\)共线,而与较大的\( \lambda\)正交的时候,\(f\)可以取最小值。
可能按照教科书,有更规范的证明方式。

发表评论

电子邮件地址不会被公开。