因为$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$可以取最小值。
可能按照教科书，有更规范的证明方式。