手撸梯度下降和python实现

一、概念与算法原理

1.1 梯度

• 从几何意义上讲，梯度就是函数变化增加最快的地方.
• 在微积分里面，对多元函数的参数求 $\partial$ 偏导数，把求得的各个参数的偏导数以向量的形式写出来，就是梯度。如函数 $f(x,y)$，其梯度向量为 $(\frac{\partial f}{\partial x},\frac{\partial f}{\partial y}{)}^T$，简称 $gradf(x,y)$ 或者 $▽f(x,y)$ 。如果是3个参数的向量梯度，就是 $(\frac{\partial f}{\partial x},\frac{\partial f}{\partial y},\frac{\partial f}{\partial z},{)}^T$ ，以此类推。

1.2 梯度下降

1.2.1 梯度下降相关概念

• 步长（Learning rate）：步长决定了在梯度下降迭代的过程中，每一步沿梯度负方向前进的长度。

• 特征（feature）：指的是样本中输入部分，比如2个样本、1个特征的数据集 { ( x ( 0 ) , y ( 0 ) ) , ( x ( 1 ) , y ( 1 ) ) } \{(x^{(0)},y^{(0)}),(x^{(1)},y^{(1)})\} {(x(0),y(0)),(x(1),y(1))}，则第一个样本特征为 $x^{(0)}$ ，第一个样本输出为 $y^{(0)}$。

• 假设函数（hypothesis function）：在监督学习中，为了拟合输入样本，而使用的假设函数，记为 $h_{\theta }(x)$ 。比如对于1个特征的n个样本 $(x^{(i)},y^{(i)})(i=1,2,...,n)$，其假设函数可写成：
  $$h_{\theta }(x)={\theta }_0+{\theta }_{1}x$$

• 损失函数（loss function）：为了评估模型拟合的好坏，通常用损失函数来度量拟合的程度。损失函数极小化，意味着拟合程度最好，对应的模型参数即为最优参数。在线性回归中，损失函数通常为样本输出和假设函数的差取平方。比如对于1个特征的n个样本 $(x^{(i)},y^{(i)})(i=1,2,...,n)$，其线性回归的损失函数为：
  $$\begin{array}{rl}J({\theta }_0,{\theta }_1)& =\sum _{i=0}^n(h_{\theta }(x{)}^{(i)}-y^{(i)}{)}^2\\ & =\sum _{i=0}^n[({\theta }_0^{(i)}+{\theta }_1^{(i)}{x}_1^{(i)})-y^{(i)}{]}^2\end{array}$$

• 梯度下降不一定能够找到全局的最优解，有可能是一个局部最优解。当然，如果损失函数是凸函数，梯度下降法得到的解就一定是全局最优解。
  

1.2.2 梯度下降的详细算法

（1）代数描述

• 1.先决条件： 确认优化模型的假设函数和损失函数。
      比如对于m个特征的线性回归，假设函数表示为 $h_{\theta }(x_1,x_2,...,x_m)={\theta }_0+{\theta }_1{x}_1+...+{\theta }_m{x}_m$，其中 ${\theta }_i(i=0,1,...,m)$ 为模型参数，$x_i(i=1,2,...,m)$为每个样本的m个特征值。为了表示可以简化，可以设$x_0=1$，假设函数则可以简化为
  $$h_{\theta }(x_0,x_1,...,x_m)=\sum _{i=0}^m{\theta }_i{x}_i$$
  对应于上面的假设函数，损失函数为：
  $$\begin{array}{r}J({\theta }_0,{\theta }_1,...,{\theta }_m)=\frac{1}{2n}\sum _{i=1}^n[h_{\theta }(x_0^{(i)},x_1^{(i)},...,x_m^{(i)})-y^{(i)}{]}^2\end{array}$$
  其中，$n$ 是样本量，$m$ 是特征数，加$\frac{1}{2n}$是为了更好地理解，也方便求偏导后抵消2次幂。

• 2.算法相关参数初始化：参数 $\theta$，步长 $\alpha$，终止距离（即梯度下降的距离）$\epsilon$

• 3.算法过程：
    1）确认梯度函数。
    2）用步长乘以损失函数的梯度，得到当前位置下降的距离。如参数 ${\theta }_i(i=0,1,...,m)$ 的梯度为 $\frac{\partial J(\theta )}{\partial {\theta }_i}$，梯度下降的距离为 $\alpha \frac{\partial J(\theta )}{\partial {\theta }_i}$，即
  $$d=\alpha \frac{1}{n}\sum _{i=1}^n[h_{\theta }(x_0^{(i)},x_1^{(i)},...,x_m^{(i)})-y^{(i)}]x^{(i)}$$
    3）确定是否所有的${\theta }_i(i=0,1,...,m)$梯度下降的距离都小于 $\epsilon$，如果小于 $\epsilon$ 则算法终止，否则进入步骤4.
    4）更新所有的${\theta }_i(i=0,1,...,m)$，更新完毕后继续转入步骤1。
  $$\begin{array}{rl}{\theta }_i& :={\theta }_i-d\\ & ={\theta }_i-\alpha \frac{1}{n}\sum _{i=1}^n[h_{\theta }(x_0^{(i)},x_1^{(i)},...,x_m^{(i)})-y^{(i)}]x^{(i)}\end{array}$$

（2）矩阵描述

• 这一部分是梯度下降法的矩阵公式，不再讲解。
• 如果需要熟悉矩阵求导建议参考张贤达的《矩阵分析与应用》一书。
  $$\begin{array}{rl}& 目标函数：h_{\theta }(X)=X\theta \\ & 损失函数：J(\theta )=\frac{1}{2}(X\theta -Y{)}^T(X\theta -Y)\\ & 损失函数偏导：\frac{\partial }{\partial \theta }J(\theta )=X^T(X\theta -Y)\\ & 梯度下降距离：\alpha X^T(X\theta -Y)\\ & 参数更新：\theta :=\theta -\alpha X^T(X\theta -Y)\end{array}$$

二、案例实战

• 设原函数 $h_{\theta }(x)={\theta }_0+{\theta }_{1}x$，我们使用以下五个数据点来拟合该函数。（正确解：$h=1+2x$）

n	1	2	3	4	5
x	1	2	3	4	5
y	3	5	7	9	11

2.1 手撸（代数形式）

  1）定义梯度公式，该案例有五个样本，两个$\theta$，所以损失函数的梯度有：
$$\begin{gathered} {\nabla }_{{\theta }_0}=\frac{1}{5}\sum _{i=1}^5[({\theta }_0^{(i)}+{\theta }_1^{(i)}{x}^{(i)})-y^{(i)}]x_0^{(i)} \\ {\nabla }_{{\theta }_1}=\frac{1}{5}\sum _{i=1}^5[({\theta }_0^{(i)}+{\theta }_1^{(i)}{x}^{(i)})-y^{(i)}]x_0^{(i)} \end{gathered}$$
其中，$x_0^{(i)}=1$.（不明白这里可以查看上文）

  2）定义梯度下降的距离公式，该案例有五个样本，两个$\theta$，所以梯度下降的距离有：
$$\begin{gathered} d0=\alpha {\nabla }_{\theta 0} \\ d1=\alpha {\nabla }_{\theta 1} \end{gathered}$$
其中，$x_0^{(i)}=1$.（不明白这里可以查看上文）

  3）定义 $\theta$ 更新公式 ：
$$\begin{gathered} {\theta }_0:={\theta }_0-d0 \\ {\theta }_1:={\theta }_1-d1 \end{gathered}$$

  4）初始化相关参数
$$\begin{array}{rl}& \alpha =0.01\\ & {\theta }_0=0,{\theta }_1=0\\ & \epsilon =0.001\end{array}$$

  5）算法过程：把$\frac{1}{n}$去掉不计算，且注意$x_0^i=1$
迭代第0次，计算出：
    两个$\theta$损失函数梯度为
$$\begin{array}{rl}{\nabla }_{{\theta }_0}& =\frac{1}{5}\ast [(0\ast 1+0\ast 2)-3]\ast 1\\ & +\frac{1}{5}\ast [(0\ast 1+0\ast 2)-5]\ast 1\\ & +\frac{1}{5}\ast [(0\ast 1+0\ast 2)-7]\ast 1\\ & +\frac{1}{5}\ast [(0\ast 1+0\ast 2)-9]\ast 1\\ & +\frac{1}{5}\ast [(0\ast 1+0\ast 2)-11]\ast 1\\ & =-0.2\ast (3+5+7+9+11)\\ & =-7\end{array}$$
$$\begin{array}{rl}{\nabla }_{{\theta }_1}& =\frac{1}{5}\ast [(0\ast 1+0\ast 2)-3]\ast 1\\ & +\frac{1}{5}\ast [(0\ast 1+0\ast 2)-5]\ast 2\\ & +\frac{1}{5}\ast [(0\ast 1+0\ast 2)-7]\ast 3\\ & +\frac{1}{5}\ast [(0\ast 1+0\ast 2)-9]\ast 4\\ & +\frac{1}{5}\ast [(0\ast 1+0\ast 2)-11]\ast 5\\ & =0.2\ast (3+10+21+46+55)\\ & =25\end{array}$$
    两个$\theta$损失函数梯度下降距离为
$$\begin{array}{rl}d0& =0.01\ast {\nabla }_{{\theta }_0}=0.01\ast (-7)=-0.07\\ d1& =0.01\ast {\nabla }_{{\theta }_1}=0.01\ast (-27)=-0.25\end{array}$$

    终止迭代判断：两个$\theta$梯度下降的距离都大于设定的阈值 $\epsilon$，继续迭代。
$$\begin{array}{rl}|d0|=0.07& >\epsilon =0.001\\ & \&\\ |d2|=0.25& >\epsilon =0.001\end{array}$$

迭代第1次，计算出：
    更新两个$\theta$
$$\begin{gathered} {\theta }_0:={\theta }_0-d0=0+0.07=0.07 \\ {\theta }_1:={\theta }_1-d1=0+0.25=0.25 \end{gathered}$$

    两个$\theta$损失函数梯度为
$$\begin{array}{rl}{\nabla }_{{\theta }_0}& =\frac{1}{5}\ast [(0.07\ast 1+0.25\ast 2)-3]\ast 1\\ & +\frac{1}{5}\ast [(0.07\ast 1+0.25\ast 2)-5]\ast 1\\ & +\frac{1}{5}\ast [(0.07\ast 1+0.25\ast 2)-7]\ast 1\\ & +\frac{1}{5}\ast [(0.07\ast 1+0.25\ast 2)-9]\ast 1\\ & +\frac{1}{5}\ast [(0.07\ast 1+0.25\ast 2)-11]\ast 1\\ & =-0.2\ast (2.423+4.43+6.43+8.43+10.43)\\ & =-6.4286\end{array}$$
$$\begin{array}{rl}{\nabla }_{{\theta }_0}& =\frac{1}{5}\ast [(0.07\ast 1+0.25\ast 2)-3]\ast 1\\ & +\frac{1}{5}\ast [(0.07\ast 1+0.25\ast 2)-5]\ast 2\\ & +\frac{1}{5}\ast [(0.07\ast 1+0.25\ast 2)-7]\ast 3\\ & +\frac{1}{5}\ast [(0.07\ast 1+0.25\ast 2)-9]\ast 4\\ & +\frac{1}{5}\ast [(0.07\ast 1+0.25\ast 2)-11]\ast 5\\ & =-0.2\ast (2.423+8.86+19.29+33.72+52.15)\\ & =-23.289\end{array}$$

    两个$\theta$损失函数梯度下降距离为
$$\begin{array}{rl}d0& =0.01\ast {\nabla }_{{\theta }_0}=0.01\ast (-6.4286)=-0.064286\\ d1& =0.01\ast {\nabla }_{{\theta }_1}=0.01\ast (-23.289)=-0.23289\end{array}$$

    终止迭代判断：两个$\theta$梯度下降的距离都大于设定的阈值 $\epsilon$，继续迭代。
$$\begin{array}{rl}|d0|=0.064286& >\epsilon =0.001\\ & \&\\ |d2|=0.23289& >\epsilon =0.001\end{array}$$

迭代第2次，计算出：
    更新两个$\theta$
$$\begin{gathered} {\theta }_0:={\theta }_0-d0=0.07+0.064286=0.134286 \\ {\theta }_1:={\theta }_1-d1=0.25+0.23289=0.48289 \end{gathered}$$

......暂且撸到这里😄

2.2 python 实现

2.2.1 代数形式

import numpy as np

def GetH_Y(x,y,Theta):
    '''
        功能：获取损失函数梯度表达式的部分：h-y
        传入：X        --> 特征，如单特征2个样本为 [[1, 1], [1, 2]]，其中x0==1
             Y        --> 对应的样本输出
             Theta    --> 模型参数
        返回:
            表达式部分：h-y 
        
    '''
    m = len(x) # m个特征，包含x0==1这个特征
    res = "+".join(['Theta[{}]*x[{}]'.format(i,i) for i in range(m)])
    
    return eval('({}) - y'.format(res))

def GetDistanceAndTheta(X,Y,Alpha,Theta):
    '''
        功能：获取梯度下降的距离和迭代更新后的theta
        传入：X       --> 特征，如单特征2个样本为 [[1, 1], [1, 2]]，其中x0==1
             Y       --> 对应的样本输出
             Theta   --> 模型参数
             Alpha   --> 步长
        返回:
            distance --> 梯度下降的距离
            Theta    --> 更新后的theta
    '''
    m = len(X[0]) # m个特征，包含x0==1这个特征
    n = len(X)    # n个样本
    distance = [] # 梯度下降距离
    for i in range(m):
        nabla = 1/n * np.sum([GetH_Y(x,y,Theta)*x[i] for x,y in zip(X,Y)]) # theta对应的损失函数梯度
        distance.append(nabla*Alpha)                                       # 梯度下降距离
    
    print("   theta: ",Theta)    # 打印目前的theta
    print("distance: ",distance) # 打印梯度下降距离
    print("---------")
    
    Theta = list(map(lambda x: x[0]-x[1], zip(Theta, distance))) # 更新 theta
    
    return distance,Theta
    
    
def GradientDescent(X,Y,Theta,Alpha,Epsilon,Freq):
    '''
        功能：使用梯度下降法计算线性函数的参数近似值（代数法）
        传入：X        --> 特征，如单特征2个样本为 [[1, 1], [1, 2]]，其中x0==1
             Y        --> 对应的样本输出
             Theta    --> 模型参数
             Alpha    --> 步长
             Freq     --> 迭代次数
             Epsilon  --> 终止迭代的阈值
        返回:
        
    '''
    print("迭代第%s次"%Freq)
    distance,Theta = GetDistanceAndTheta(X,Y,Alpha,Theta) # 开始迭代
    AbsDistance = list(map(lambda x:abs(x),distance))     # 取梯度下降距离的绝对值
    
    # 判断任何一个theta的梯度下降距离是否大于给定的阈值，如果为True，则继续迭代，否则终止迭代.
    while (np.array(AbsDistance) > Epsilon).any():
        Freq += 1                                             # 迭代次数加 1
        print("迭代第%s次"%Freq)                               # 打印迭代次数
        distance,Theta = GetDistanceAndTheta(X,Y,Alpha,Theta) # 继续迭代
        AbsDistance = list(map(lambda x:abs(x),distance))     # 取梯度下降距离的绝对值
    
    print("迭代终止!")
    
if __name__ == "__main__":
    x = [[1,1], [1,2], [1,3], [1,4], [1,5]] # 其中x0==1
    y = [3, 5, 7, 9, 11]
    
    freq    = 0      # 迭代次数初始化为0
    theta   = [0,0]  # 参数初始化
    alpha   = 0.01   # 步长
    epsilon = 0.0001 # 迭代终止条件的阈值
    GradientDescent(x,y,theta,alpha,epsilon,freq)    

$........................................................................$

  不难看出，前几次迭代的结果与上文手撸版一致（结果可能存在小数点差异）；最终的$\theta$值也很接近最优解${\theta }_0=1,{\theta }_1=2$。
  上文计算$\theta$的损失函数梯度时加了$\frac{1}{n}$，其实加和不加对对结果影响不大。比如将程序中的代码去掉$\frac{1}{n}$

改为

运行结果

  最终迭代了419次（加$\frac{1}{n}$迭代了1151次），只是迭代次数减少了（针对本案例而言），但结果是相差不大的。

2.2.2 矩阵形式

import numpy as np

def GetDistanceAndTheta(X,Y,Alpha,Theta):
    '''
        功能：获取梯度下降的距离和迭代更新后的theta
        传入：X       --> 特征，如单特征2个样本为 [[1, 1], [1, 2]]，其中x0==1
             Y       --> 对应的样本输出
             Theta   --> 模型参数
             Alpha   --> 步长
        返回:
             D       --> 梯度下降的距离
             Theta   --> 更新后的theta
    '''
    D = np.dot(Alpha*np.array(X).T,np.dot(X,Theta)-Y)

    print("   Theta: ",Theta)  # 打印目前的theta
    print("Distance: ",D)      # 打印梯度下降距离
    print("---------")
    
    Theta -= D # 更新 theta
    
    return D,Theta
    
    
def GradientDescent(X,Y,Theta,Alpha,Epsilon,Freq):
    '''
        功能：使用梯度下降法计算线性函数的参数近似值（代数法）
        传入：X        --> 特征，如单特征2个样本为 [[1, 1], [1, 2]]，其中x0==1
             Y        --> 对应的样本输出
             Theta    --> 模型参数
             Alpha    --> 步长
             Freq     --> 迭代次数
             Epsilon  --> 终止迭代的阈值
        返回:
        
    '''
    print("迭代第%s次"%Freq)
    D,Theta = GetDistanceAndTheta(X,Y,Alpha,Theta) # 迭代第0次
    AbsD = list(map(lambda x:abs(x),D))            # 取梯度下降距离的绝对值
    
    # 判断任何一个theta的梯度下降距离是否大于给定的阈值，如果为True，则继续迭代，否则终止迭代.
    while (np.array(AbsD) > Epsilon).any():
        Freq += 1                                      # 迭代次数加 1
        print("迭代第%s次"%Freq)                         # 打印迭代次数
        D,Theta = GetDistanceAndTheta(X,Y,Alpha,Theta) # 继续迭代
        AbsD = list(map(lambda x:abs(x),D))            # 取梯度下降距离的绝对值
        AbsD = list(map(lambda x:abs(x),D))            # 取梯度下降距离的绝对值
        
    print("迭代终止!")
    
if __name__ == "__main__":
    x = [[1,1], [1,2], [1,3], [1,4], [1,5]] # 其中x0==1
    y = [3, 5, 7, 9, 11]
    
    freq    = 0      # 迭代次数初始化为0
    theta   = [0,0]  # 参数初始化
    alpha   = 0.01   # 步长
    epsilon = 0.0001 # 迭代终止条件的阈值
    GradientDescent(x,y,theta,alpha,epsilon,freq)    

$........................................................................$

• 该结果与代数形式下不加 $\frac{1}{n}$ 的结果一致。
• 矩阵形式的写法代码较为简洁、运行速度较快.

参考：梯度下降（Gradient Descent）小结