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

假设我们要拟合一个函数，目前我们知道的函数值是 ( x 1 , y 1 ) , ( x 2 , y 2 ) , ⋯   , ( x n , y n ) (x_1,y_1),(x_2,y_2),cdots,(x_n,y_n) (x1​,y1​),(x2​,y2​),⋯,(xn​,yn​)，而我们要求 y = f ( x ) y=f(x) y=f(x)，使得它的函数值与真实值的偏差 r i = f ( x i ) − y i r_i=f(x_i)-y_i ri​=f(xi​)−yi​的平方和 ∑ i = 1 n r i 2 = ∑ i = 1 n ( f ( x i ) − y i ) 2 sumlimits_{i=1}^nr_i^2=sumlimits_{i=1}^n(f(x_i)-y_i)^2 i=1∑n​ri2​=i=1∑n​(f(xi​)−yi​)2最小。

假如我们已经选定了一类函数（模型）去拟合，比如 f ( x ) = e a x + b f(x)=e^{ax+b} f(x)=eax+b，现在要确定其中参数的值。设参数为 a 1 , a 2 , ⋯   , a m a_1,a_2,cdots,a_m a1​,a2​,⋯,am​，将 f f f表示为 f ( x , a 1 , a 2 , ⋯   , a m ) f(x,a_1,a_2,cdots,a_m) f(x,a1​,a2​,⋯,am​)。损失函数，即偏差的平方和为 Q ( a 1 , a 2 , ⋯   , a m ) = ∑ i = 1 n ( f ( x i , a 1 , a 2 , ⋯   , a m ) − y i ) 2 Q(a_1,a_2,cdots,a_m)=sumlimits_{i=1}^n(f(x_i,a_1,a_2,cdots,a_m)-y_i)^2 Q(a1​,a2​,⋯,am​)=i=1∑n​(f(xi​,a1​,a2​,⋯,am​)−yi​)2。让 Q Q Q取得极值，就需要让 Q Q Q对 a 1 , a 2 , ⋯   , a m a_1,a_2,cdots,a_m a1​,a2​,⋯,am​的偏导数都为 0 0 0。对于参数 a k a_k ak​， Q Q Q对它的偏导数为 ∂ Q ∂ a k = 2 ∑ i = 1 n [ ( f ( x i , a 1 , a 2 , ⋯   , a m ) − y i ) ∂ f ∂ a k ∣ x i ] frac{partial Q}{partial a_k}=2sumlimits_{i=1}^nleft[left(f(x_i,a_1,a_2,cdots,a_m)-y_iright)left.frac{partial f}{partial a_k}right|_{x_i}right] ∂ak​∂Q​=2i=1∑n​[(f(xi​,a1​,a2​,⋯,am​)−yi​)∂ak​∂f​∣∣∣∣​xi​​]所以就是要使得 ∀ k = 1 , 2 , ⋯   , m , ∑ i = 1 n [ ( f ( x i , a 1 , a 2 , ⋯   , a m ) − y i ) ∂ f ∂ a k ∣ x i ] = 0 forall k=1,2,cdots,m,\sumlimits_{i=1}^nleft[left(f(x_i,a_1,a_2,cdots,a_m)-y_iright)left.frac{partial f}{partial a_k}right|_{x_i}right]=0 ∀k=1,2,⋯,m,i=1∑n​[(f(xi​,a1​,a2​,⋯,am​)−yi​)∂ak​∂f​∣∣∣∣​xi​​]=0

下面以一次函数拟合为例。设 f ( x ) = a x + b f(x)=ax+b f(x)=ax+b，则 Q ( a , b ) = ∑ i = 1 n ( a x i + b − y i ) 2 Q(a,b)=sumlimits_{i=1}^n(ax_i+b-y_i)^2 Q(a,b)=i=1∑n​(axi​+b−yi​)2， ∂ Q ∂ a = 2 ∑ i = 1 n ( a x i + b − y i ) x i = 0 ① ∂ Q ∂ b = 2 ∑ i = 1 n ( a x i + b − y i ) = 0 ② begin{aligned}frac{partial Q}{partial a}&=2sumlimits_{i=1}^n(ax_i+b-y_i)x_i=0qquad&①\frac{partial Q}{partial b}&=2sumlimits_{i=1}^n(ax_i+b-y_i)=0&②end{aligned} ∂a∂Q​∂b∂Q​​=2i=1∑n​(axi​+b−yi​)xi​=0=2i=1∑n​(axi​+b−yi​)=0​①②​由 ② ② ②得 a ∑ i = 1 n x i + n b = ∑ i = 1 n y i asumlimits_{i=1}^nx_i+nb=sumlimits_{i=1}^ny_i ai=1∑n​xi​+nb=i=1∑n​yi​令 x ‾ = ∑ i = 1 n x i n overline x=frac{sumlimits_{i=1}^n x_i}{n} x=ni=1∑n​xi​​， y ‾ = ∑ i = 1 n y i n overline y=frac{sumlimits_{i=1}^n y_i}{n} y​=ni=1∑n​yi​​，则有 a x ‾ + b = y ‾ ③ aoverline x+b=overline yqquad ③ ax+b=y​③由 ① ① ①得 a ∑ i = 1 n x i 2 + b ∑ i = 1 n x i = ∑ i = 1 n x i y i asumlimits_{i=1}^n x_i^2+bsumlimits_{i=1}^n x_i=sumlimits_{i=1}^n x_iy_i ai=1∑n​xi2​+bi=1∑n​xi​=i=1∑n​xi​yi​令 x 2 ‾ = ∑ i = 1 n x i 2 n , x y ‾ = ∑ i = 1 n x i y i n overline{x^2}=frac{sumlimits_{i=1}^n x_i^2}{n},overline{xy}=frac{sumlimits_{i=1}^n x_iy_i}{n} x2=ni=1∑n​xi2​​,xy​=ni=1∑n​xi​yi​​，则有 a x 2 ‾ + b x ‾ = x y ‾ ④ aoverline{x^2}+boverline{x}=overline{xy}qquad④ ax2+bx=xy​④ ③ × x ‾ ③timesoverline x ③×x得 a x ‾ 2 + b x ‾ = x ˉ y ˉ ⑤ a{overline x}^2+boverline x=bar xbar yqquad⑤ ax2+bx=xˉyˉ​⑤ ⑤ − ④ ⑤-④ ⑤−④得 a ( x ‾ 2 − x 2 ‾ ) = x ˉ y ˉ − x y ‾ a({overline x}^2-overline{x^2})=bar xbar y-overline{xy} a(x2−x2)=xˉyˉ​−xy​即 a = x ˉ y ˉ − x y ‾ x ‾ 2 − x 2 ‾ a=frac{bar xbar y-overline{xy}}{{overline x}^2-overline{x^2}} a=x2−x2xˉyˉ​−xy​​由 ③ ③ ③得 b = y ‾ − a x ‾ b=overline y-aoverline x b=y​−ax

一次函数拟合的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()

执行结果：

时间复杂度： Θ ( n ) Theta(n) Θ(n)。