【高等数学笔记】彻底弄懂最小二乘法（Least Squares Method）

假设我们要拟合一个一元函数，目前我们知道的自变量及因变量的值是$(x_1,y_1),(x_2,y_2),\cdots ,(x_n,y_n)$，而我们要求$y=f(x)$，使得用它作为近似函数时，函数值与真实值的偏差$r_i=f(x_i)-y_i$的平方和$\sum _{i=1}^n{r}_i^2=\sum _{i=1}^n(f(x_i)-y_i{)}^2$最小。

假如我们已经选定了一类函数（模型）去拟合，现在要确定其中参数的值。比如我们选择了$f(x)=e^{ax+b}$这种模型（一般用于拟合人口增长），就要确定参数$a,b$的值。设参数为$a_1,a_2,\cdots ,a_m$，将$f$表示为$f(x,a_1,a_2,\cdots ,a_m)$。损失函数，即偏差的平方和为$Q(a_1,a_2,\cdots ,a_m)=\sum _{i=1}^n(f(x_i,a_1,a_2,\cdots ,a_m)-y_i{)}^2$。注意，这里$Q$是参数$a_1,a_2,\cdots ,a_m$的函数，而不是$x,y$的函数。让$Q$取得极值，就需要让$Q$对$a_1,a_2,\cdots ,a_m$的偏导数都为$0$。对于参数$a_k$，$Q$对它的偏导数为$$\frac{\partial Q}{\partial a_k}=2\sum _{i=1}^n\left[\left(f(x_i,a_1,a_2,\cdots ,a_m)-y_i\right){\frac{\partial f}{\partial a_k}|}_{x_i}\right]$$让它等于$0$，就是要使得$$\begin{gathered} \forall k=1,2,\cdots ,m, \\ \sum _{i=1}^n\left[\left(f(x_i,a_1,a_2,\cdots ,a_m)-y_i\right){\frac{\partial f}{\partial a_k}|}_{x_i}\right]=0 \end{gathered}$$

---

下面以一次函数拟合为例。设$f(x)=ax+b$，要通过已知的数据点$(x_1,y_1),(x_2,y_2),\cdots ,(x_n,y_n)$来确定参数$a,b$的值。损失函数$Q(a,b)=\sum _{i=1}^n(a{x}_i+b-y_i{)}^2$，$$\begin{array}{rlr}\frac{\partial Q}{\partial a}& =2\sum _{i=1}^n(a{x}_i+b-y_i)x_i=0& ①\\ \frac{\partial Q}{\partial b}& =2\sum _{i=1}^n(a{x}_i+b-y_i)=0& ②\end{array}$$由$②$得$$a\sum _{i=1}^n{x}_i+nb=\sum _{i=1}^n{y}_i$$令$\overline{x}=\frac{\sum _{i=1}^n{x}_i}{n}$，$\overline{y}=\frac{\sum _{i=1}^n{y}_i}{n}$（分别为$x,y$的平均值），则有$$a\overline{x}+b=\overline{y}③$$由$①$得$$a\sum _{i=1}^n{x}_i^2+b\sum _{i=1}^n{x}_i=\sum _{i=1}^n{x}_i{y}_i$$令$\overline{x^2}=\frac{\sum _{i=1}^n{x}_i^2}{n},\overline{xy}=\frac{\sum _{i=1}^n{x}_i{y}_i}{n}$，则有$$a\overline{x^2}+b\overline{x}=\overline{xy}④$$$③\times \overline{x}$得$$a{\overline{x}}^2+b\overline{x}=\bar{x}\bar{y}⑤$$$⑤-④$得$$a({\overline{x}}^2-\overline{x^2})=\bar{x}\bar{y}-\overline{xy}$$由此算出直线的斜率$$a=\frac{\bar{x}\bar{y}-\overline{xy}}{{\overline{x}}^2-\overline{x^2}}$$由$③$得$$b=\overline{y}-a\overline{x}$$

---

一次函数拟合的Python代码实现：

# Least Squares Method (Linear)
# Author: seh_sjij

import matplotlib.pyplot as plt
import numpy as np

class LeastSquareMethod(object):
    def __init__(self, x, y):
        self.x = x
        self.y = y
        self.n = len(self.x)
        if len(self.y) != self.n:
            raise Exception(
                'LeastSquareMethod: len(x) != len(y)')

    def Calculate(self):
        self.xmean = self.xsquaremean = 0
        # average of x and x^2
        for x_i in self.x:
            self.xmean += x_i
            self.xsquaremean += x_i * x_i
        self.xmean /= self.n
        self.xsquaremean /= self.n
        self.ymean = 0
        # average of y
        for y_i in self.y:
            self.ymean += y_i
        self.ymean /= self.n
        self.xymean = 0
        # average of xy
        for i in range(0, self.n):
            self.xymean += self.x[i] * self.y[i]
        self.xymean /= self.n
        a = (self.xmean * self.ymean - self.xymean) \
            / (self.xmean * self.xmean - self.xsquaremean)
        b = self.ymean - a * self.xmean
        return (a, b) # y = ax + b

if __name__ == '__main__':
    x = [1, 2, 4, 6, 7, 9, 10, 12]
    y = [0.7, 2.25, 4.64, 5.69, 7.40, 8.57, 10.72, 11.64]
    lsm = LeastSquareMethod(x, y)
    a, b = lsm.Calculate()
    plt.scatter(x, y)
    plt.plot(x, [a * x_i + b for x_i in x])
    plt.title('Least Squares Method: y = %fx + %f' % (a, b))
    plt.show()

执行结果：

时间复杂度：$\Theta (n)$。