引言

众所周知的是，在大学课程中一般只会教授一种拟合方法(也即参数估计方法)——最小二乘法。这是一种直接求解的方法，非常的有效，不仅是损失最小解，而且是最大似然解。只不过，有一个缺点，它只能解决线性方程参数问题，对于非线性曲线，就无能为力了。大部分情况下还是将其转换成线性问题，再使用最小二乘法。

然而，并非所有的问题都能转换为线性问题，甚至并非所有目标建模公式的参数都能有解析解，其他学科如机器学习等学科如何解决这个参数估计问题？答案是——《最优化方法》，其中最常用的是梯度下降法，不去寻找解析解，而是寻找其导数/梯度。因为导数/梯度具有如下优点

1. 导数/梯度永远指向数值变动最快的方向（梯度的性质）
2. 导数/梯度（除非有边界和断点）永远指向一个到达极值的路径

我们的主要目标是
arg min ⁡ θ L ( θ ) = 1 n ∑ i = 1 n ( y i − f θ ( x i ) ) 2 \argmin_\theta L(\theta)=\dfrac{1}{n}\sum_{i=1}^{n} (y_i-f_\theta(x_i))^2 θargmin​L(θ)=n1​i=1∑n​(yi​−fθ​(xi​))2
这是最小二乘问题，也叫均方误差最小化问题。其中 θ \theta θ就是我们想要求的参数。 f θ ( x ) f_\theta(x) fθ​(x)就是我们的模型，即带有参数 θ \theta θ的函数。 x i , y i x_i,y_i xi​,yi​分别是数据中第 i i i对自变量和因变量。 L ( θ ) L(\theta) L(θ)是损失函数，用来衡量我们拟合的效果。
而梯度下降法的主要内容是
θ ← θ − η ∇ θ L u n t i l ∇ θ L = 0 \theta \gets \theta- \eta\nabla_\theta L \quad until \quad \nabla_\theta L=0 θ←θ−η∇θ​Luntil∇θ​L=0
其中 η \eta η是学习率，目的是不要让 θ \theta θ变动幅度过大（特别是一些有限制的量，比如log函数的底数），导致不能收敛。 ∇ θ \nabla_\theta ∇θ​是求对 θ \theta θ的梯度。取负号是因为梯度默认是增大函数值的方向，这里是最小化问题，所以取其相反数。由这里可知，梯度大小不是我们所需要的，我们只要梯度的方向。

问题引入

关于logistics增长问题，是这样的，假设某生物种群数量记为 N N N，最开始为常数 N 0 N_0 N0​。
假设有未知参数 r r r， K K K，其中 r r r是固定增长率，表明在没有任何限制的情况下种群数量的增长速度。 K K K是环境容纳量，是表明某个环境对生物的限制导致的最大容纳数量。

它们的关系如下
d N d t = r × N × ( 1 − N K ) \dfrac{dN}{dt}=r\times N\times(1-\dfrac{N}{K}) dtdN​=r×N×(1−KN​)
我们的数据是按照时间 t i t_i ti​记录的种群数量 N i N_i Ni​。我们的目标就是求出 r r r， K K K。
这是一条微分方程。关于它的原函数可以参考我写的另一篇文章
它的图像是这样的
其中 N 0 = 10 , r = 1 , K = 1000 N_0=10,r=1,K=1000 N0​=10,r=1,K=1000。横坐标是时间 t t t，纵坐标是种群数量 N N N

链式求导法则

首先，我们先明确实际输出与目标输出，按照logistic公式，我们的输出
f i = f ( t i , N i ) = r × N i × ( 1 − N i K ) f_i=f(t_i,N_i)=r\times N_i\times(1-\dfrac{N_i}{K}) fi​=f(ti​,Ni​)=r×Ni​×(1−KNi​​)
我们的目标
d N ( t i ) d t i ≈ N i + 1 − N i t i + 1 − t i \dfrac{dN(t_i)}{dt_i}\approx \dfrac{N_{i+1}-N_i}{t_{i+1}-t_i} dti​dN(ti​)​≈ti+1​−ti​Ni+1​−Ni​​
如果我们更细究一点， d N ( t ξ i ) d t ξ i = N i + 1 − N i t i + 1 − t i , t ξ i ∈ [ t i , t i + 1 ] \dfrac{dN(t_{\xi_i})}{dt_{\xi_i}}= \dfrac{N_{i+1}-N_i}{t_{i+1}-t_i},t_{\xi_i}\in[t_i,t_{i+1}] dtξi​​dN(tξi​​)​=ti+1​−ti​Ni+1​−Ni​​,tξi​​∈[ti​,ti+1​]。我们可以大约令 N ξ i = N i + 1 + N i 2 , t ξ i = t i + 1 + t i 2 N_{\xi_i}=\frac{N_{i+1}+N_i}{2},t_{\xi_i}=\frac{t_{i+1}+t_i}{2} Nξi​​=2Ni+1​+Ni​​,tξi​​=2ti+1​+ti​​，将点 ( t ξ i , N ξ i ) (t_{\xi_i},N_{\xi_i}) (tξi​​,Nξi​​)对应导数认为 N ξ i ′ ≈ N i + 1 − N i t i + 1 − t i N_{\xi_i}'\approx\dfrac{N_{i+1}-N_i}{t_{i+1}-t_i} Nξi​′​≈ti+1​−ti​Ni+1​−Ni​​。
那么损失函数
L ( r , K ) = 1 n − 1 ∑ i = 1 n − 1 ( N ξ i ′ − f ( t ξ i , N ξ i ) ) 2 = 1 n − 1 ∑ i = 1 n − 1 ( N ξ i ′ − f ξ i ) 2 L(r,K)=\dfrac{1}{n-1}\sum_{i=1}^{n-1} (N_{\xi_i}'-f(t_{\xi_i},N_{\xi_i}))^2=\dfrac{1}{n-1}\sum_{i=1}^{n-1} (N_{\xi_i}'-f_{\xi_i})^2 L(r,K)=n−11​i=1∑n−1​(Nξi​′​−f(tξi​​,Nξi​​))2=n−11​i=1∑n−1​(Nξi​′​−fξi​​)2

在logistics增长问题中，未知量是 r r r， K K K，所以我们的目标是求出 ∇ r L , ∇ K L \nabla_r L,\nabla_K L ∇r​L,∇K​L。
∇ r L = ∂ L ∂ r , ∇ K L = ∂ L ∂ K \nabla_r L=\dfrac{\partial L}{\partial r},\quad\nabla_K L=\dfrac{\partial L}{\partial K} ∇r​L=∂r∂L​,∇K​L=∂K∂L​
根据链式求导法则，我们有

∂ L ∂ r = ∑ i = 1 n − 1 ∂ L ∂ f ξ i ∂ f ξ i ∂ r \dfrac{\partial L}{\partial r}=\sum_{i=1}^{n-1}\dfrac{\partial L}{\partial f_{\xi_i}}\dfrac{\partial f_{\xi_i}}{\partial r} ∂r∂L​=i=1∑n−1​∂fξi​​∂L​∂r∂fξi​​​
同理 K K K也有类似形式。
而我们有
∂ L ∂ f ξ i = − 2 ( N ξ i ′ − f ξ i ) , ∂ f ξ i ∂ r = N ξ i ( 1 − N ξ i K ) , ∂ f ξ i ∂ K = r N ξ i 2 K 2 \dfrac{\partial L}{\partial f_{\xi_i}}=-2(N_{\xi_i}'-f_{\xi_i}),\quad\dfrac{\partial f_{\xi_i}}{\partial r}=N_{\xi_i}(1-\dfrac{N_{\xi_i}}{K}),\quad\dfrac{\partial f_{\xi_i}}{\partial K}=r\dfrac{N_{\xi_i}^2}{K^2} ∂fξi​​∂L​=−2(Nξi​′​−fξi​​),∂r∂fξi​​​=Nξi​​(1−KNξi​​​),∂K∂fξi​​​=rK2Nξi​2​​
所以
∇ r L = ∂ L ∂ r = − 2 ∑ i = 1 n − 1 ( N ξ i ′ − f ( t ξ i , N ξ i ) ) N ξ i ( 1 − N ξ i K ) ∇ K L = ∂ L ∂ K = − 2 ∑ i = 1 n − 1 ( N ξ i ′ − f ( t ξ i , N ξ i ) ) r N ξ i 2 K 2 \nabla_r L=\dfrac{\partial L}{\partial r}=-2\sum_{i=1}^{n-1}(N_{\xi_i}'-f(t_{\xi_i},N_{\xi_i}))N_{\xi_i}(1-\dfrac{N_{\xi_i}}{K})\\ \nabla_K L=\dfrac{\partial L}{\partial K}=-2\sum_{i=1}^{n-1}(N_{\xi_i}'-f(t_{\xi_i},N_{\xi_i}))r\dfrac{N_{\xi_i}^2}{K^2}\\ ∇r​L=∂r∂L​=−2i=1∑n−1​(Nξi​′​−f(tξi​​,Nξi​​))Nξi​​(1−KNξi​​​)∇K​L=∂K∂L​=−2i=1∑n−1​(Nξi​′​−f(tξi​​,Nξi​​))rK2Nξi​2​​

实践

首先是数据的生成

import numpy as np
N0=10
K_real=1000
r_real=1
c=math.log(N0/(K_real-N0))#原函数
def ground_true(t):return(K_real/(1+np.exp(-r_real*t-c)))x=np.arange(0,10,0.1)
y=ground_true(x)# dN/dt 目标输出
y_diff = (y[1:]-y[:-1])/(x[1:]-x[:-1])# N_xi 处理后的输入
y_xi = (y[1:]+y[:-1])/2# 这里不想用n-1，直接以n表示了
n = y_xi.shape[0]

原数据效果就是一开始给的图，而经过处理的数据则是如下效果

import matplotlib.pyplot as plt
plt.plot(x,y)
plt.plot(x[:-1],y_xi)
plt.show()

定义必要的函数

# 模型
def f(r, K, N):return r*N*(1-N/K)# loss函数    
def loss_of_rK(r, K):return 1/n*np.sum((y_diff-f(r, K, y_xi))**2)# ∂L/∂r
def diff_r(r, K):return -2.0/n*np.sum((y_diff-f(r, K, y_xi))*y_xi*(1-y_xi/K))# ∂L/∂K
def diff_K(r, K):return -2.0/n*np.sum((y_diff-f(r, K, y_xi))*r*(y_xi**2)/(K**2))

在选择梯度下降法的起点时，尽量靠近最终目标，收敛时间才短。
前面说过，梯度的大小不是我们所关心的内容，所以我加入了一个归一化函数 arctan ⁡ \arctan arctan（归一化函数还有很多，如sigmoid等等）用来取出它的方向，用参数本身的绝对值当做步伐大小。（其实这比较像是强行把L2损失改为L1损失）

#初始化
rx = np.random.rand()
Kx = np.random.rand()*np.max(y_xi)#选接近最大值为初始值比较合理# 学习率
eta = 1e-1
# 误差精度
epsilon= 1
while (loss_of_rK(rx, Kx) > epsilon):rx = rx-eta*np.arctan(diff_r(rx, Kx))*np.abs(rx)Kx = Kx-eta*np.arctan(diff_K(rx, Kx))*np.abs(Kx)

迭代了41610步，结果为

它的拟合轨迹如下，横轴为r，纵轴为K

需要注意的是，梯度下降法没有办法保证解是最优/解存在/解可达。也没有办法保证收敛。它只是一种凸优化手段。尽管如此，机器学习等学科仍然大量的使用它，俗称炼丹。

本文链接：https://my.lmcjl.com/post/2887.html