简介
简单介绍了最小二乘法原理,以直线拟合为例结合原理出发推导了拟合算法过程,最后c++实现拟合算法。
最小二乘法
最小二乘法(英文:least square method)是一种常用的数学优化方法,所谓二乘就是平方的意思。这平方一词指的是在拟合一个函数的时候,通过最小化误差的平方来确定最佳的匹配函数,所以最小二乘、最小平方指的就是拟合的误差平方达到最小。
推导过程
以直线拟合为例,已知有一组平面上的点集$(x_1, y_1), (x_2,y_2),...(x_n,y_n)$。基于这些点拟合一条直线,设直线方程为
$$y = ax + b$$
那么那么算法的输入就是这些点集,需要求取的为直线方程的参数a,b。
(1)平方偏差之和为
$$S_{\epsilon^2} = \sum_{i = 1}^{n}{(y_i - y)}^2$$
$$ = \sum_{i = 1}^{n}{(y_i - (ax_i + b))}^2$$
那么,以上公式已知的是$x_i, y_i$, 未知且要求取的是a、b。不同的a、b会得到不同的$S_{\epsilon^2}$,求取的是在$S_{\epsilon^2}$最小的时候求取a、b。这是一个二元(a,b)函数,此问题实际上就是多元函数的极值与最值问题,要求解函数极值时二元变量数值,这里要用到二元函数取极值的必要条件
$$设z = f(x,y)在点(x_0, y_0) \begin{cases} \text{一阶偏导数存在} \\ \text{取极值}\end{cases}$$
则
$$f_x^,(x_0,y_0)=0$$
以及
$$f_y^,(x_0,y_0)=0$$
那么对 $S_{\epsilon^2}$ 求偏导且使得偏导为0,此时 $S_{\epsilon^2}$取得极值点最小值
$$\frac{\partial}{\partial{a}}{S_{\epsilon^2}} = \sum_{i=1}^n{2(y_i -(ax_i + b)(-x_i))=0}$$
$$ = \sum_{i=1}^n{(ax_i^2 + bx_i - y_ix_i)=0}$$
$$ => (\sum_{i=1}^n{x_i^2})a + (\sum_{i=1}^n x_i)b = \sum_{i=1}^n{y_ix_i}$$
$$\frac{\partial}{\partial{b}}S_{\epsilon^2} = \sum_{i=1}^n2(y_i -(ax_i + b)(-1) = 0$$
$$= \sum_{i=1}^n(ax_i + b - y_i) = 0$$
$$ => (\sum_{i=1}^nx_i)a + nb = \sum_{i =1}^n y_i$$
联立两式可得方程组
$$\begin{cases}(\sum_{i=1}^n{x_i^2})a + (\sum_{i=1}^n x_i)b = \sum_{i=1}^n{y_ix_i} \\ (\sum_{i=1}^nx_i)a + nb = \sum_{i =1}^n y_i\end{cases}$$
求解方程组可得a、b的值。将方程组转为矩阵形式表示
$${\begin{pmatrix} \sum_{i=1}^n{x_i^2} & \sum_{i=1}^n x_i \\ \sum_{i=1}^nx_i & n \end{pmatrix}} {\begin{pmatrix} a \\ b \end{pmatrix}} = \begin{pmatrix}{\sum_{i=1}^ny_ix_i} \\ {\sum_{i=1}^ny_i} \end{pmatrix}$$
倘若给点集的每个点加上权重,方程组则变为
$$\begin{cases}(\sum_{i=1}^n{x_i^2}w_i)a + (\sum_{i=1}^n x_iw_i)b = \sum_{i=1}^n{y_ix_iw_i} \\ (\sum_{i=1}^nx_iw_i)a + {\sum_{i=1}^n w_i} = \sum_{i =1}^n y_iw_i\end{cases}$$
矩阵表示
$${\begin{pmatrix} \sum_{i=1}^n{x_i^2w_i} & \sum_{i=1}^n x_iw_i \\ \sum_{i=1}^nx_iw_i & {\sum_{i=1}^n w_i} \end{pmatrix}} {\begin{pmatrix} a \\ b \end{pmatrix}} = \begin{pmatrix}{\sum_{i=1}^ny_ix_iw_i} \\ {\sum_{i=1}^ny_iw_i} \end{pmatrix}$$
实现代码
实现思路
(1)输入点集$(x_1, y_1), (x_2,y_2),...(x_n,y_n)$
(2)根据点集数据和权重系数分别计算出
$${\begin{pmatrix} \sum_{i=1}^n{x_i^2w_i} & \sum_{i=1}^n x_iw_i \\ \sum_{i=1}^nx_iw_i & {\sum_{i=1}^n w_i} \end{pmatrix}}$$
以及
$$ \begin{pmatrix}{\sum_{i=1}^ny_ix_iw_i} \\ {\sum_{i=1}^ny_iw_i} \end{pmatrix}$$
矩阵。
(3)使用OpenCV解矩阵函数cv::solve()
计算出
$${\begin{pmatrix} a \\ b \end{pmatrix}}$$
下面放出函数实现,完整代码已托管在github,如果觉得有对你有所帮助,不妨给个star吧。
代码链接:https://github.com/mangosroom/machine-vision-algorithms-library/tree/master/src/linefit
void mv::LineFit::FitLineByRegression()
{
// 设置权重 | weights setting
_weigths = std::vector<double>(_points.size(), 1);
// AX = B
// 构造A矩阵 | Construct A mat
const int N = 2;
cv::Mat A = cv::Mat::zeros(N, N, CV_64FC1);
for (int row = 0;row < A.rows;row++)
{
for (int col = 0; col < A.cols;col++)
{
for (int k = 0;k < _points.size();k++)
{
A.at<double>(row, col) = A.at<double>(row, col) + pow(_points[k].x, row + col) * _weigths[k];
}
}
}
//构造B矩阵 | Construct B mat
cv::Mat B = cv::Mat::zeros(N, 1, CV_64FC1);
for (int row = 0;row < B.rows;row++)
{
for (int k = 0;k < _points.size();k++)
{
B.at<double>(row, 0) = B.at<double>(row, 0) + pow(_points[k].x, row)*_points[k].y * _weigths[k];
}
}
// 求解A*X = B | Solve the A*X = B
cv::Mat X;
cv::solve(A, B, X,cv::DECOMP_LU);
// y = b + ax
_result.b = X.at<double>(0,0);
_result.a = X.at<double>(1, 0);
}//FitLineByRegression
Reference
本文由芒果浩明发布,转载请注明出处。
本文链接:https://mangoroom.cn/opencv/least-square-method-line-fit.html