opencv最小二乘拟合求圆心 最小二乘拟合 python

目录

• 1 最小二乘拟合（方法一）

• 1.1 数学推导
• 1.2 算例
• 1.3 Python 代码

• 2.最小二乘拟合（方法二）

• 2.1 数学推导
• 2.2 算例
• 2.3 Python 代码

• 3 最小二乘法拟合（方法三）

• 3.1 数学推导
• 3.2 算例
• 3.3 Python 代码

• 4 利用sklearn.linear_model()

• 4.1 参考资料
• 4.2 Python 代码

1 最小二乘拟合（方法一）

本章介绍的是我在上研究生课程《数值分析》时所学的内容，也是最一般的解法，可以对任意函数进行拟合。

1.1 数学推导

对一组数据$(x_i,y_i),i=1,2\cdots,m$，要在某个函数类$\Phi=\{\phi_0(x),\phi_1(x),\cdots,\phi_n(x)\},n\ll m$中构造一个函数$\phi(x)=\sum_{k=0}^{n}a_k\phi_k(x)$，使得$I=\sum_{i=1}^{m}\left[ \phi(x_i)-y_i\right]^2$取得极小值。
此处的函数类$\Phi$就是指很多种类函数的集合，例如幂函数、三角函数、指数函数、有理函数、多项式函数等。例如，常用的多项式函数类有$\Phi_1=\{1,x,x^2,\cdots,x^n\}$，$\Phi_2=\{1,\underbrace{x,x,\cdots,x}_{n}\}$。其中$\Phi_1$为n阶多项式，$\Phi_2$为n元一阶多项式，分别对应着多项式拟合和多维线性回归。
残差函数$I$定义为
$I(a_0,a_1,\cdots,a_n)=\sum_{i=1}^{m}\left[a_0\phi_0(x_i)+a_1\phi_1(x_i)+\cdots+a_n\phi_n(x_i) -y_i\right]^2$
令$I$对系数$a_k$的偏导为0，即可求出使残差函数$I$最小的系数$a_k$的值
$\begin{aligned} \frac{\partial I}{\partial a_k}&=2\sum_{i=1}^m\left[ \phi(x_i)-y_i \right]\frac{\partial \phi(x_i)}{\partial a_k} \\ &=2\sum_{i=1}^m \left[ 2\sum_{j=0}^n a_j\phi_j(x_i) - y_i\right]\phi_k(x_i)\\ &=2\left\{ \sum_{j=0}^n a_j\left[\sum_{i=1}^m \phi_j(x_i)\phi_k(x_i)\right] - \sum_{i=1}^m y_i\phi_k(x_i) \right\} \end{aligned}$
简便起见，记
$\begin{aligned} (\phi_j,\phi_k)&=\sum_{i=1}^m \phi_j(x_i)\phi_k(x_i)\\ (y,\phi_k)&=\sum_{i=1}^m y_i\phi_k(x_i) \end{aligned}$
则$\frac{\partial I}{\partial a_k}$可写成如下形式
$\frac{\partial I}{\partial a_k} = 2\left[\sum_{j=0}^na_j(\phi_j,\phi_k) - (y,\phi_k)\right] = 0$
化简为
$\sum_{j=0}^na_j(\phi_j,\phi_k) = (y,\phi_k)$
注意到，公式中的系数$a_k$中$k=0,1,\cdots,n$，则实际上的操作步骤为分别对每个系数求偏导，并令偏导为0，从而求出使残差函数$I$取到极小值的所有的系数。写成矩阵的形式
$\begin{bmatrix} (\phi_0,\phi_0)& (\phi_0,\phi_1)& \cdots & (\phi_0,\phi_n) \\ (\phi_1,\phi_0)& (\phi_1,\phi_1)&\cdots&(\phi_1,\phi_n)\\ \vdots &\vdots&&\vdots\\ (\phi_n,\phi_0)& (\phi_n,\phi_1)&\cdots&(\phi_n,\phi_n) \end{bmatrix} \begin{bmatrix} a_0\\a_1\\ \vdots \\ a_n \end{bmatrix}= \begin{bmatrix} (y,\phi_0)\\(y,\phi_1)\\ \vdots \\(y,\phi_n) \end{bmatrix}$
上式方程组称为法方程组（Normal Equation），向量$\begin{bmatrix}a_0\\a_1\\ \vdots \\ a_n\end{bmatrix}$称为回归系数（Regression Coefficients）。用更一般的形式将上式再次改写
$Ax=b$
则可通过矩阵求逆的方式求出回归系数
$x = A^{-1}b$

1.2 算例

利用一阶多项式$\Phi_1=\{1,x\}$通过1.1节的方法对如下数据进行最小二乘拟合。

Independent variable x	Dependent variable y
0	394.33
4	329.50
8	291.00
12	255.17
16	229.33
20	204.83
24	179.00
28	163.83
32	150.33

则法方程组为

$\begin{bmatrix} 9&\sum x_i\\\sum x_i&\sum x_i^2 \end{bmatrix} \begin{bmatrix} a_0\\a_1 \end{bmatrix}= \begin{bmatrix} \sum y_i \\ \sum x_iy_i \end{bmatrix}$

带入数据

$\begin{bmatrix} 9&144\\144&3264 \end{bmatrix} \begin{bmatrix} a_0\\a_1 \end{bmatrix}= \begin{bmatrix} 2197.32 \\ 28167.72 \end{bmatrix}$

求解得到

$\begin{bmatrix} a_0\\a_1 \end{bmatrix}= \begin{bmatrix} 360.636667\\-7.280625 \end{bmatrix}$

则可得到线性回归方程

$y = -7.280625x + 360.636667$

拟合结果如下

1.3 Python 代码

#%% import neccesary packages

import numpy as np
import matplotlib.pyplot as plt
#%% generate the dataset

x = np.arange(0,33,4)
y = np.array([394.33, 329.50, 291.00, 255.17, 229.33, 204.83, 179.00, 163.83, 150.33])
print(x)
#%% calculate the normal equation

A11 = 9
A12 = np.sum(x)
A21 = A12
A22 = np.sum(np.power(x,2))
A = np.array([[A11,A12],[A21,A22]])
print(A)

b1 = np.sum(y)
b2 = np.sum(x*y)
b=np.array([b1, b2])
print(b)
#%% solve the polynamial coefficient a0 & a1

a = np.linalg.inv(A).dot(b)
a0 = a[0]
a1 = a[1]
print(a)
#%% plot the picture
fig, ax = plt.subplots(1, 1, figsize=(6, 3))
ax.scatter(x, y, color='black')
plt.xlabel("x")
plt.ylabel("y")
plt.grid()
ax.plot(x, a0 + a1 * x)
plt.show()

---

2.最小二乘拟合（方法二）

本章介绍我自己的推导思路，本文以多项式拟合为例说明思路，其他函数的拟合思路相同。

2.1 数学推导

对一组数据$(x_i,y_i),i=1,2\cdots,m$进行多项式拟合，假设多项式为
$f(x)=a_nx^n+a_{n-1}x^{n-1}+\cdots+a_1x+a_0$
对任意数据点$(x_k,y_k)$，它与多项式$f(x)$在$y$轴方向上的距离称为残差$\delta_k$，定义为
$\delta_k=f(x_k)-y_k$
拟合的目标是寻找一组多项式系数$\{a_0,a_1,\cdots,a_n\}$使得残差函数$I$取最小值，本章将直接利用最小二乘法求解。将每个数据点的残差写成矩阵形式
$\begin{bmatrix} 1& x_1& x^2_1&\cdots & x^n_1 \\ 1&x_2& x^2_2&\cdots&x^n_2\\ \vdots &\vdots &\vdots&&\vdots\\ 1&x_m& x^2_m&\cdots&x^n_m \end{bmatrix} \begin{bmatrix} a_0\\a_1\\ \vdots \\ a_n \end{bmatrix}- \begin{bmatrix} y_1\\y_2\\ \vdots \\y_m \end{bmatrix} = \begin{bmatrix} \delta_1\\\delta_2\\ \vdots \\\delta_m \end{bmatrix}$
定义残差函数$I(a_0,a_1,\cdots,a_n)$为
$I=\parallel\delta_i\parallel ^2_2=\parallel f(x_i) - y_i \parallel ^2_2$
要使残差函数$I$最小，本质上是求如下方程的最小二乘解
$\begin{bmatrix} 1& x_1& x^2_1&\cdots & x^n_1 \\ 1&x_2& x^2_2&\cdots&x^n_2\\ \vdots &\vdots &\vdots&&\vdots\\ 1&x_m& x^2_m&\cdots&x^n_m \end{bmatrix} \begin{bmatrix} a_0\\a_1\\ \vdots \\ a_n \end{bmatrix}= \begin{bmatrix} y_1\\y_2\\ \vdots \\y_m \end{bmatrix}$
上式可写成一般形式
$Ax=b$
（1）当$A$为方阵时，其最小二乘解为
$x=-(A^TA)^{-1}A^Tb$
（2）在大多数情况下，$A$不为方阵，可用Moore-Penrose逆$A^+$求解
定理1： 设$A\in \mathbb{C} ^{m \times n}$是秩为$r(r>0)$的矩阵，且$A$的满秩分解为
$A=FG\quad (F\in \mathbb{C} ^{m \times r},G\in \mathbb{C} ^{r \times n}秩都为r)$
则
$A^+=G^H(GG^H)^{-1}(F^HF)^{-1}F^H$
定理2： 设$A\in \mathbb{C} ^{m \times n},b\in \mathbb{C} ^{m}$有解的充要条件是
$AA^+b=b$
其通解为
$x=A^+b+(I-A^+A)y$
其唯一极小范数最小二乘解为
$x=A^+b$

2.2 算例

本章算例与1.2节相同，以作对照。
利用一阶多项式$y=a_0+a_1x$对如下数据进行最小二乘拟合。

Independent variable x	Dependent variable y
0	394.33
4	329.50
8	291.00
12	255.17
16	229.33
20	204.83
24	179.00
28	163.83
32	150.33

列写$Ax=b$的一般形式

$\begin{bmatrix} 1& 0\\ 1& 4\\ 1& 8\\ 1& 12\\ 1& 16\\ 1& 20\\ 1& 24\\ 1& 28\\ 1& 32 \end{bmatrix} \begin{bmatrix} a_0\\ a_1\\ \end{bmatrix}= \begin{bmatrix} 394.33\\ 329.50\\ 291.00\\ 255.17\\ 229.33\\ 204.83\\ 179.00\\ 163.83\\ 150.33 \end{bmatrix}$

利用$Python$命令$numpy.linalg.pinv()$求解矩阵$A$的$A^+$逆，并求其唯一极小范数最小二乘解

$\begin{bmatrix} a_0\\a_1 \end{bmatrix}= \begin{bmatrix} 360.636667\\-7.280625 \end{bmatrix}$

求解结果与第一章相同。

2.3 Python 代码

#%% import neccesary packages

import numpy as np
import matplotlib.pyplot as plt
#%% generate the dataset

x = np.arange(0,33,4)
y = np.array([394.33, 329.50, 291.00, 255.17, 229.33, 204.83, 179.00, 163.83, 150.33])
print(x)
#%% solve the polynamial coefficient a0 & a1

A = np.vstack((np.array([1]*9), x))
print(A)
pinv_A = np.linalg.pinv(A).T
print(pinv_A)
a = pinv_A.dot(y)
print(a)
a0 = a[0]
a1 = a[1]
#%% plot the picture

fig, ax = plt.subplots(1, 1, figsize=(6, 3))
ax.scatter(x, y, color='black')
plt.xlabel("x")
plt.ylabel("y")
plt.grid()
ax.plot(x, a0 + a1 * x)
plt.show()

---

3 最小二乘法拟合（方法三）

本章重点介绍线性回归算法。其实该算法本质上和第一章方法一样，只不过将其限定为一阶多项式，推导出更直观的表达方法。

3.1 数学推导

假设线性回归模型如下
$y_i = \beta_0+\beta_1x_i+\epsilon_i$
式中$\epsilon_i$为模型噪声服从高斯分布：$\epsilon \sim N(0,\sigma^2)$。
假设需要对一组数据$(x_i,y_i),i=1,2\cdots,m$进行线性回归，定义残差函数
$Q = \sum_{i=1}^n \left[ y_i-(\beta_0+\beta_1x_i) \right]^2$
优化目标为最小化残差函数
$\mathop{min}\limits_{\beta_0,\beta_1}: Q$
和第一章方法一样，令残差函数$Q$对系数$\beta_0、\beta_1$的偏导为0
$\begin{aligned} \frac{\partial Q}{\partial \beta_0}=-2\sum_{i=1}^{n}\left[ y_i-(\beta_0+\beta_1x_i) \right]&=0 \\ \Rightarrow \qquad n\beta_0+\beta_1\sum_{i=1}^nx_i &=\sum_{i=1}^ny_i \qquad (1) \end{aligned}$
$\begin{aligned} \frac{\partial Q}{\partial \beta_1}=-2\sum_{i=1}^{n}x_i\left[ y_i-(\beta_0+\beta_1x_i) \right]&=0 \\ \Rightarrow \qquad \beta_0\sum_{i=1}^nx_i+\beta_1\sum_{i=1}^nx_i^2 &=\beta_1\sum_{i=1}^nx_iy_i \qquad (2) \end{aligned}$
联立式（1）和式（2）解方程得到$\beta_0、\beta_1$的代数表达式
$\hat{\beta}_0=\frac{(\sum_{i=1}^nx_i^2)(\sum_{i=1}^ny_i)-(\sum_{i=1}^nx_i )(\sum_{i=1}^nx_iy_i)}{ n\sum_{i=1}^nx_i^2-(\sum_{i=1}^nx_i)^2}$
$\hat{\beta}_1=\frac{\sum_{i=1}^nx_iy_i-(\sum_{i=1}^ny_i)(\sum_{i=1}^nx_i )}{ n\sum_{i=1}^nx_i^2-(\sum_{i=1}^nx_i)^2}$
为表达简便，记
$S_{xy}=\sum_{i=1}^nx_iy_i-\frac{1}{n}(\sum_{i=1}^nx_i)(\sum_{i=1}^ny_i)=\sum_{i=1}^n(x_i-\overline{x})(y_i-\overline{y})$
$S_{xx}=\sum_{i=1}^nx_i^2-\frac{1}{n}(\sum_{i=1}^nx_i)^2=\sum_{i=1}^n(x_i-\overline{x})^2$
$S_{yy}=\sum_{i=1}^ny_i^2-\frac{1}{n}(\sum_{i=1}^ny_i)^2=\sum_{i=1}^n(y_i-\overline{y})^2$
其中，$\overline{x}、\overline{y}$分别为拟合数据$x_i、y_i$的平均值，$i=1,2,\cdots,m$
$\begin{aligned} \overline{x}&= \frac{1}{n}\sum_{i=1}^nx_i \\ \overline{y}&= \frac{1}{n}\sum_{i=1}^ny_i \end{aligned}$
则$\beta_0、\beta_1$的表达式可改写为
$\begin{aligned} \hat{\beta_0}&= \overline{y}-\beta_1\overline{x} \\ \hat{\beta_1}&= \frac{S_{xy}}{S_{xx}} \end{aligned}$

3.2 算例

本章算例与1.2节相同，以作对照。
对如下数据进行线性回归。

Independent variable x	Dependent variable y
0	394.33
4	329.50
8	291.00
12	255.17
16	229.33
20	204.83
24	179.00
28	163.83
32	150.33

计算系数$\hat{\beta}_0、\hat{\beta}_1$

$\begin{aligned} \overline{x}&=16 \\ \overline{y}&=244.15 \\ S_{xy}&=-6989.40 \\ S_{xx}&=960 \\ \hat{\beta_1}&= \frac{S_{xy}}{S_{xx}} = -7.280625 \\ \hat{\beta_0}&= \overline{y}-\beta_1\overline{x}=360.636667 \end{aligned}$

计算结果和第一章、第二章结果一致。

$y = 360.636667-7.280625x$

3.3 Python 代码

#%% import neccesary packages

import numpy as np
import matplotlib.pyplot as plt
#%% generate the dataset

x = np.arange(0,33,4)
y = np.array([394.33, 329.50, 291.00, 255.17, 229.33, 204.83, 179.00, 163.83, 150.33])
print(x)
#%%

x_bar = np.sum(x)/9
y_bar = np.sum(y)/9
Sxy = np.sum((x-x_bar)*(y-y_bar))
Sxx = np.sum(np.power(x-x_bar,2))
beta_1 = Sxy/Sxx
beta_0 = y_bar - x_bar * beta_1
print(beta_0)
print(beta_1)
#%%

fig, ax = plt.subplots(1, 1, figsize=(6, 3))
ax.scatter(x, y, color='black')
plt.xlabel("x")
plt.ylabel("y")
plt.grid()
ax.plot(x, beta_0 + beta_1 * x)
plt.show()

---

4 利用sklearn.linear_model()

4.1 参考资料

官网链接：Scikit-learn官网链接 教学1：python-sklearn学习笔记 第一节 linear_model 教学2：numpy中newaxis函数的基本用法及其理解（傻瓜版）

4.2 Python 代码

#%% import neccesary packages

import numpy as np
import matplotlib.pyplot as plt
import sklearn.linear_model as lm

#%% generate the dataset

x = np.arange(0,33,4)
y = np.array([394.33, 329.50, 291.00, 255.17, 229.33, 204.83, 179.00, 163.83, 150.33])
print(x)
#%%

x_lr = np.linspace(0,40,500)
lr = lm.LinearRegression()
lr.fit(x[:, np.newaxis],y)
y_lr = lr.predict(x_lr[:, np.newaxis])
#%%

fig, ax = plt.subplots(1, 1, figsize=(6, 3))
ax.scatter(x, y, color='black')
plt.xlabel("x")
plt.ylabel("y")
plt.grid()
ax.plot(x_lr, y_lr)