Python实现传统梯度下降法（一）

就算中二也要好好学习
2,142 阅读10分钟

1.梯度下降法简介

梯度下降法是一种求函数极小值以及其极小值点的算法，是一种迭代优化的算法，所谓迭代优化，就是根据当前函数值的情况，根据一定的规则，调整寻优的方向和步长，最终求得极小值。 梯度下降法认为，在对函数进行寻优的过程中，函数梯度的反方向是函数下降最快的方向，就像是在爬山过程中，想迅速下山（不考虑人的身体状况）则应该向最陡峭的方向行进，负梯度就是这个最陡峭的方位。

梯度下降相关公式

设y=f(x)y=f(x)为目标函数，grady=∇f(x)grady=\nabla f(x)，每次下降的步长为η\eta。 假设我们给定的迭代初值为x0x_0，则它的下一个迭代值x1x_1为： x1=x0−η∇f(x)x0x_1=x_0-\eta \nabla f(x)_{x_0}

同理可知，不失一般性，xl=xl−1−η∇f(x)xl−1x_{l}=x_{l-1}-\eta \nabla f(x)_{x_{l-1}},其中l为迭代次数。 当梯度下降到低于一个阈值或者函数值差值下降至低于一个阈值时，停止迭代。

以上是梯度下降法的简要介绍，以下是对梯度下降法的Python实现。同时我们讨论梯度下降的相关参数对梯度下降效果的影响。

#导入相关库
import matplotlib.pyplot as plt
import numpy as np

2.简单的一元的梯度下降函数

首先我们编辑了一个简单的医院的梯度下降函数，其主要参数为初始值x0x_0,以及学习率η\eta。我们将此函数命名为gradient1()。我们设置初始精度为10−510^{-5}，迭代停止条件为梯度绝对值（范数）低于精度值。

最终输出为：

x1:求解值

n_iter:迭代次数

finaly:最终函数取值

finaldy:最终梯度取值

lf:迭代过程中函数的取值

lx:迭代过程中x的取值

def gradient1(x0,eta,func,dfunc,eps=1e-5):
    '''
    args:
    
    x0:初值
    eta:学习率
    eps:精度,默认为1e-5
    func:需要求的函数
    dfunc:导函数 
    
    returns:
    x1:最终迭代值
    count:迭代次数
    func(x1):终值
    dfunc(x1):最终梯度值
    fval:迭代值列表
    xval:参数值列表
    '''
    count=0
    epsilon=abs(dfunc(x0))
    fval=[func(x0)]
    xval=[x0]
    x1=x0
    while epsilon>eps:
        x1=x0-eta*dfunc(x0)
        x0=x1
        epsilon=abs(dfunc(x0))
        fval.append(func(x1))
        xval.append(x1)
        count+=1
    return {'x1':x1,'n_iter':count,'finaly':func(x1),'finaldy':dfunc(x1),'lf':fval,'lx':xval}

2.一元梯度下降简单结果呈现

2.1凸函数

2.1.1一元二次函数

我们设定一个一元二次函数y=x2+4x+6y=x^2+4x+6，其导函数为dy=2x+4dy=2x+4， 通过求导计算我们可以知道，该函数的极小值点为-2，极小值为2。以下是该函数的设定函数以及几何图像。

def qf1(x):
    '''
    定义一个一元二次函数
    x:输入的参数
    '''
    return x**2+4*x+6
def dqf1(x):
    '''
    定义以上一元二次函数的导数形式
    x:输入的参数
    '''
    return 2*x+4

xsam=np.linspace(-4,0,100)
ysam=[qf1(i) for i in xsam]
plt.plot(xsam,ysam)
plt.scatter(-2,qf1(-2),s=100,color='red')
plt.show()

output_7_0.png

2.1.2学习率

我们首先观察学习率对于梯度下降迭代效果的影响，我们设置了一组学习率，分别使用这些学习率进行迭代，得到的结果如下：

diceta=dict()
for eta1 in [0.05,0.1,0.2,0.4,0.5,0.7,0.85,0.95]:
    diceta[eta1]=gradient1(0,eta1,qf1,dqf1,eps=1e-5)

首先，我们可以看到，随着学习率的增加，迭代次数不是越多越好也不是越少越好，它存在一个先下降后上升的趋势。

对于任何一个学习率下的梯度下降，得到的极小值点的结果是相近的，都在精度之内。

for eta1 in [0.05,0.1,0.2,0.4,0.5,0.7,0.85,0.95]:
    print('当学习率为{}时，梯度下降法迭代次数为{}'.format(eta1,diceta[eta1]['n_iter']))    

当学习率为0.05时，梯度下降法迭代次数为123
当学习率为0.1时，梯度下降法迭代次数为58
当学习率为0.2时，梯度下降法迭代次数为26
当学习率为0.4时，梯度下降法迭代次数为9
当学习率为0.5时，梯度下降法迭代次数为1
当学习率为0.7时，梯度下降法迭代次数为15
当学习率为0.85时，梯度下降法迭代次数为37
当学习率为0.95时，梯度下降法迭代次数为123

for eta1 in [0.05,0.1,0.2,0.4,0.5,0.7,0.85,0.95]:
    print('当学习率为{}时，梯度下降法解为{}'.format(eta1,diceta[eta1]['x1']))    

当学习率为0.05时，梯度下降法解为-1.9999952917593056
当学习率为0.1时，梯度下降法解为-1.9999952109514347
当学习率为0.2时，梯度下降法解为-1.9999965883654365
当学习率为0.4时，梯度下降法解为-1.999998976
当学习率为0.5时，梯度下降法解为-2.0
当学习率为0.7时，梯度下降法解为-2.000002147483648
当学习率为0.85时，梯度下降法解为-2.000003712423184
当学习率为0.95时，梯度下降法解为-2.0000047082406947

plt.rcParams['font.sans-serif'] = ['fangsong'] # 步骤一（替换sans-serif字体）
plt.rcParams['font.size'] = 12 #设置字体大小
plt.rcParams['axes.unicode_minus'] = False   # 步骤二（解决坐标轴负数的负号显示问题）
fig,ax=plt.subplots(2,4,figsize=(16,8))
etaa=[0.05,0.1,0.2,0.4,0.5,0.7,0.85,0.95]
cmap1=plt.get_cmap('Set1_r')
colors=cmap1(np.arange(8))
for j in range(2):
    for i in range(4):
        eta1=etaa[j*4+i]
        ax[j][i].plot(xsam,ysam)
        ax[j][i].scatter(-2,qf1(-2),s=100,color='red')
        ax[j][i].plot(diceta[eta1]['lx'],diceta[eta1]['lf'],c=colors[j*4+i])
        label='初值为0$\eta$={}时的迭代情况'.format(eta1)
        ax[j][i].set_title(label)
        ax[j][i].text(x=-3,y=4,s='迭代次数{}'.format(diceta[eta1]['n_iter']))
plt.show()

从下图我们可以看到，当学习率较小时，迭代一步一步沿着同一个方向进行，当梯度增大到一定程度的时候，我们可以看到，迭代值在极值左右发生了震荡，且学习率越大，震荡次数越多。

output_12_0.png

2.1.3初值

我们使用不同的初值，同一个学习率η=0.2\eta=0.2得到结果如下：

可以初步看出，在本部分函数的迭代中，初值对于函数迭代的影响并不大。

dicx0=dict()
for x0 in [-3,-2.5,-1,0]:
    dicx0[x0]=gradient1(x0,0.2,qf1,dqf1,eps=1e-5)

for x0 in [-3,-2.5,-1,0]:
    print('当初值为{}时，梯度下降法迭代次数为{}'.format(x0,dicx0[x0]['n_iter']))   

当初值为-3时，梯度下降法迭代次数为24
当初值为-2.5时，梯度下降法迭代次数为23
当初值为-1时，梯度下降法迭代次数为24
当初值为0时，梯度下降法迭代次数为26

for x0 in [-3,-2.5,-1,0]:
    print('当初值为{}时，梯度下降法解为{}'.format(x0,dicx0[x0]['x1']))   

当初值为-3时，梯度下降法解为-2.0000047383813384
当初值为-2.5时，梯度下降法解为-2.0000039486511154
当初值为-1时，梯度下降法解为-1.9999952616186618
当初值为0时，梯度下降法解为-1.9999965883654365

2.2其他函数

以上我们使用的函数为严格凸的函数，它有且只有一个极小值点，因此最终总能得到一个结果，但是如果我们需要迭代的函数并不是一个严格凸的函数呢，如果其存在多个局部极小值点呢，此处我们定义的函数如下：

f(x)=x4−4x2−xf(x)=x^4-4x^2-x

df(x)=4x3−8x−1df(x)=4x^3-8x-1

def sif1(x):
    '''
    定义一个一元二次函数
    x:输入的参数
    '''
    return x**4-4*x**2-x
def dsif1(x):
    '''
    定义以上一元二次函数的导数形式
    x:输入的参数
    '''
    return 4*x**3-8*x-1

xsam4=np.linspace(-2,2,100)
ysam4=[sif1(i) for i in xsam4]
plt.plot(xsam4,ysam4)
plt.show()

以下是其几何图像： output_19_0.png

2.2.1初值为2

首先我们使用初值为2，并改换学习率，得到的结果如下：

diceta2=dict()
for eta1 in [0.01,0.05,0.1,0.2]:
    diceta2[eta1]=gradient1(2,eta1,sif1,dsif1,eps=1e-3)

for eta1 in [0.01,0.05,0.1,0.2]:
    print('当学习率为{}时，梯度下降法迭代次数为{}'.format(eta1,diceta2[eta1]['n_iter']))

在运行过程中，我发现，当学习率超过0.2（未必精确）的时候，迭代次数过多，程序停止。

当学习率为0.01时，梯度下降法迭代次数为44
当学习率为0.05时，梯度下降法迭代次数为6
当学习率为0.1时，梯度下降法迭代次数为25
当学习率为0.2时，梯度下降法迭代次数为14007

for eta1 in [0.01,0.05,0.1,0.2]:
    print('当学习率为{}时，结果为{}'.format(eta1,diceta2[eta1]['x1']))

还可以看到，在学习率为0.01,0.05,0.1时，迭代的最终结果都是1.47,然而当学习率为0.2时，迭代的最终结果为-1.34，猜测由于学习率较大，导致迭代飞跃到另一局部极小值的位置出不来了（突发恶疾）。

当学习率为0.01时，结果为1.4730464134091468
当学习率为0.05时，结果为1.47298744179518
当学习率为0.1时，结果为1.472949438559799
当学习率为0.2时，结果为-1.3469689471855626

fig1,ax11=plt.subplots(1,4,figsize=(16,4))
for i in range(4):
    eta1=[0.01,0.05,0.1,0.2]
    ax11[i].plot(xsam4,ysam4)
    ax11[i].plot(diceta2[eta1[i]]['lx'],diceta2[eta1[i]]['lf'])
    ax11[i].set_title('初值为2$\eta$={}时的迭代情况'.format(eta1[i]))

从下图，可以看出，当学习率为0.2的时候，迭代确实出现了这样的问题。

output_24_0.png

2.2.2初值为-2

当初值为-2的时候，得到的结论类似，具体可以看以下的迭代次数，和迭代图像。

diceta22=dict()
for eta1 in [0.01,0.05,0.1,0.2]:
    diceta22[eta1]=gradient1(-2,eta1,sif1,dsif1,eps=1e-3)

for eta1 in [0.01,0.05,0.1,0.2]:
    print('当学习率为{}时，梯度下降法迭代次数为{}'.format(eta1,diceta22[eta1]['n_iter']))    

当学习率为0.01时，梯度下降法迭代次数为57
当学习率为0.05时，梯度下降法迭代次数为9
当学习率为0.1时，梯度下降法迭代次数为12
当学习率为0.2时，梯度下降法迭代次数为12853

for eta1 in [0.01,0.05,0.1,0.2]:
    print('当学习率为{}时，结果为{}'.format(eta1,diceta22[eta1]['x1']))

当学习率为0.01时，结果为-1.3470670587886961
当学习率为0.05时，结果为-1.346961120419268
当学习率为0.1时，结果为-1.3470589888762234
当学习率为0.2时，结果为-1.3470495983375566

fig1,ax11=plt.subplots(1,4,figsize=(16,4))
for i in range(4):
    eta1=[0.01,0.05,0.1,0.2]
    ax11[i].plot(xsam4,ysam4)
    ax11[i].plot(diceta22[eta1[i]]['lx'],diceta22[eta1[i]]['lf'])
    ax11[i].set_title('初值为2$\eta$={}时的迭代情况'.format(eta1[i]))

output_29_0.png

3.推广到n维空间（以二维为例）

3.1改写梯度下降形式

和一元的梯度下降函数结构基本相同，但是在比较梯度精度的时候，将梯度的绝对值该车梯度的l2范数。

def gradient2(x0,eta,func,dfunc,eps=1e-5):
    '''
    args:
    
    x0:初值，numpy数组形式
    eta:学习率
    eps:精度,默认为1e-5
    func:需要求的函数
    dfunc:导函数 
    
    returns:
    x1:最终迭代值
    count:迭代次数
    func(x1):终值
    dfunc(x1):最终梯度值
    fval:迭代值列表
    xval:参数值列表
    '''
    count=0
    epsilon=sum(abs(dfunc(x0))**2)**0.5
    fval=[func(x0)]
    xval=[x0]
    x1=x0
    while epsilon>eps:
        x1=x0-eta*dfunc(x0)
        x0=x1
        epsilon=sum(abs(dfunc(x0))**2)**0.5
        fval.append(func(x1))
        xval.append(x1)
        count+=1
    return {'x1':x1,'n_iter':count,'finaly':func(x1),'finaldy':dfunc(x1),'lf':fval,'lx':xval}

3.2定义二元函数

我们定义了函数z=x2+y2z=x^2+y^2

其梯度为[2x,2y][2x,2y]

以下是这个函数的等高线图和3D图像。

#定义一个椭圆函数
def eryuan(x):
    '''
    x是numpy数组
    '''
    return x[0]**2+x[1]**2
def deryuan(x):
    d1=2*x[0]
    d2=2*x[1]
    return np.array([d1,d2])

from matplotlib import cm
from matplotlib.ticker import LinearLocator

xsam3 = np.arange(-20, 20, 0.25)
ysam3 = np.arange(-20, 20, 0.25)
X, Y = np.meshgrid(xsam3, ysam3)
Z = X**2+Y**2
levels = np.linspace(np.min(Z), np.max(Z), 20)
fig, ax = plt.subplots()
ax.contour(X, Y, Z, levels=levels,cmap=cm.coolwarm)
plt.show()

output_36_0.png

fig, ax = plt.subplots(subplot_kw={"projection": "3d"},figsize=(10,10))
X = np.arange(-5, 5, 0.25)
Y = np.arange(-5, 5, 0.25)
X, Y = np.meshgrid(X, Y)
Z = X**2 + Y**2
surf = ax.plot_surface(X, Y, Z, cmap=cm.PiYG,
                       linewidth=0, antialiased=False)
ax.set_zlim(0,25)
ax.zaxis.set_major_locator(LinearLocator(10))
ax.zaxis.set_major_formatter('{x:.02f}')
fig.colorbar(surf)
plt.show()

output_37_0.png

3.3结果展示

多元梯度下降的结果展示如下

resulter=gradient2(np.array([10,10]),0.4,eryuan,deryuan,eps=1e-5)

resulter

{'x1': array([1.024e-06, 1.024e-06]),
 'n_iter': 10,
 'finaly': 2.097151999999989e-12,
 'finaldy': array([2.048e-06, 2.048e-06]),
 'lf': [200,
  8.0,
  0.31999999999999984,
  0.012799999999999987,
  0.0005119999999999991,
  2.0479999999999953e-05,
  8.191999999999976e-07,
  3.2767999999999874e-08,
  1.3107199999999944e-09,
  5.242879999999976e-11,
  2.097151999999989e-12],
 'lx': [array([10, 10]),
  array([2., 2.]),
  array([0.4, 0.4]),
  array([0.08, 0.08]),
  array([0.016, 0.016]),
  array([0.0032, 0.0032]),
  array([0.00064, 0.00064]),
  array([0.000128, 0.000128]),
  array([2.56e-05, 2.56e-05]),
  array([5.12e-06, 5.12e-06]),
  array([1.024e-06, 1.024e-06])]}

xsam3 = np.arange(-20, 20, 0.25)
ysam3 = np.arange(-20, 20, 0.25)
X, Y = np.meshgrid(xsam3, ysam3)
Z = X**2+Y**2
levels = np.linspace(np.min(Z), np.max(Z), 20)
fig, ax = plt.subplots(figsize=(9,9))
ax.contour(X, Y, Z, levels=levels,cmap=cm.coolwarm)
ax.plot(np.array(resulter['lx']).T[0],np.array(resulter['lx']).T[1],marker='*',c='red')
plt.show()

output_40_0.png

avatar
文章
阅读
粉丝