1、多项式曲线拟合之最小二乘法

1.1 问题来源

1801年，意大利天文学家朱赛普·皮亚齐发现了第一颗小行星谷神星。经过40天的跟踪观测后，由于谷神星运行至太阳背后，使得皮亚齐失去了谷神星的位置。随后全世界的科学家利用皮亚齐的已有观测数据开始寻找谷神星，但是根据大多数人计算的结果来寻找谷神星都没有结果。只有时年24岁的高斯所计算的谷神星的轨道，被奥地利天文学家海因里希·奥尔伯斯的观测所证实，使天文界从此可以预测到谷神星的精确位置。同样的方法也产生了哈雷彗星等很多天文学成果。高斯使用的方法就是最小二乘法，该方法发表于1809年他的著作《天体运动论》中。

1.2 数学本质

采用最小二乘法进行曲线拟合的本质是通过样本集构造范德蒙德矩阵，将一元n次多项式非线性回归问题转化为n元一次线性回归问题。

给定一组数据点 p_i(x_i,y_i)$ ,其中 i=1,2,...m 。求近似曲线 $y=\varphi(x)$ ,使其与 y=f(x) 的偏差最小。

常见的曲线拟合方法：

使偏差绝对值之和最小

$\mathop{min}_{\varphi}\sum_{i=1}^m{\left|\delta_i\right|} = \sum_{i=1}^m{\left|\varphi(x_i)-y_i\right|}$

使最大的偏差绝对值最小

$\mathop{min}_{\varphi}\ \mathop{max}_{i}{\left|\delta_i\right|} = \left|\varphi(x_i)-y_i\right|$

使偏差平方和最小

$\mathop{min}_{\varphi}\sum_{i=1}^m{\delta_i^2} = \sum_{i=1}^m(\varphi(x_i)-y_i)^2$

其中按照偏差平方和最小的原则选取拟合曲线，并且采取二项式方程为拟合曲线的方法，称为最小二乘法。

1.3 问题定义

$minimize \qquad \parallel{Ax-b}\parallel_2^2$

1.4 问题特性

已是一种成熟的工业技术
已有可靠和高效的算法解决此类问题
存在可解析解： $x^* = ({\mathbf{A}^\mathrm{T}A)}^{-1}{\mathbf{A}^\mathrm{T}b}$

1.5 数学推导

设定拟合多项式为： $y = a_0+a_1{x}+\dots+{a_k}{x^k}$
偏差平方和表示如下：

$R^2 = \sum_{i=1}^n[y_i - (a_0+a_1{x_i}+\dots+{a_k}{x_i}^k)]^2$

对右侧等式求 $\alpha_i$ 偏导数，以求的符合条件的 $\alpha$ 值：

对 $\alpha_0$ 求偏导：

$-2\sum_{i=1}^n[y_i - (a_0+a_1{x_i}+\dots+{a_k}{x_i}^k)] = 0$

对 $\alpha_1$ 求偏导：

$-2\sum_{i=1}^n[y_i - (a_0+a_1{x_i}+\dots+{a_k}{x_i}^k)]{x_i} = 0$

对 $\alpha_2$ 求偏导：

$-2\sum_{i=1}^n[y_i - (a_0+a_1{x_i}+\dots+{a_k}{x_i}^k)]{x_i}^2 = 0$

$\vdots$

对 $\alpha_k$ 求偏导：

$-2\sum_{i=1}^n[y_i - (a_0+a_1{x_i}+\dots+{a_k}{x_i}^k)]{x_i}^k = 0$

等式化简

$a_0{n}+a_1\sum_{i=1}^n{x_i}+\dots+a_k\sum_{i=1}^n{x_i}^k = \sum_{i=1}^n{y_i}$

$a_0\sum_{i=1}^n{x_i}+a_1\sum_{i=1}^n{x_i}^2+\dots+a_k\sum_{i=1}^n{x_i}^{k+1} = \sum_{i=1}^n{x_i}{y_i}$

$a_0\sum_{i=1}^n{x_i}^2+a_1\sum_{i=1}^n{x_i}^3+\dots+a_k\sum_{i=1}^n{x_i}^{k+2} = \sum_{i=1}^n{x_i}^2{y_i}$

$a_0\sum_{i=1}^n{x_i}^k+a_1\sum_{i=1}^n{x_i}^{k+1}+\dots+a_k\sum_{i=1}^n{x_i}^{k+k} = \sum_{i=1}^n{x_i}^k{y_i}$

矩阵表示

$\left[ \begin{array}{cccc} n & \sum_{i=1}^n{x_i} & \cdots & \sum_{i=1}^n{x_i}^k\\ \sum_{i=1}^n{x_i} & \sum_{i=1}^n{x_i}^2 & \cdots & \sum_{i=1}^n{x_i}^{k+1}\\ \vdots & \vdots & \ddots & \vdots\\ \sum_{i=1}^n{x_i}^k & \sum_{i=1}^n{x_i}^{k+1} & \cdots & \sum_{i=1}^n{x_i}^{k+k}\\ \end{array} \right] \left[ \begin{array}{cccc} a_0\\ a_1\\ \vdots\\ a_k\\ \end{array} \right] = \left[ \begin{array}{cccc} \sum_{i=1}^n{y_i}\\ \sum_{i=1}^n{x_i}{y_i}\\ \vdots\\ \sum_{i=1}^n{x_i}^k{y_i}\\ \end{array} \right]$

矩阵简化

$X= \left[ \begin{array}{cccc} 1 & x_1 & x_1^2 & \cdots x_1^k \\ 1 & x_2 & x_2^2 & \cdots x_2^k \\ \cdots & \cdots & \cdots & \cdots \\ 1 & x_n & x_n^2 & \cdots x_n^k \\ \end{array} \right] \qquad Y= \left[ \begin{array}{cccc} y_1 \\ y_2 \\ \cdots \\ y_n \\ \end{array} \right]$

上述矩阵可简化为： $\mathbf{X}^\mathrm{T}Xa=\mathbf{X}^\mathrm{T}Y$ 。