1、Nesterov动量

在动量法提出不久之后，Nesterov动量也随之被提了出来，此方法属于动量法的进一步发展，与动量法不同的是，在动量法中会用每一次迭代中计算得来动量项作为参数调整的方向和速率，而在Nesterov动量法中，在计算新的动量项之前，会用当前的动量项进行一次参数的调整，然后用调整的临时参数计算梯度值，然后再利用该梯度值进行新的动量项的计算，最后进行参数的更新，此过程可以看作在动量法中添加了一个较正因子。Nesterov动量法的核心思想在于在每次计算梯度的时候使用下一次的权重，然后将下一次的变化情况反映在这次计算的梯度上。
Nesterov动量法的具体步骤如下：

使用当前的动量项进行参数的临时调整；
使用调整后的临时参数计算出临时梯度值g；
使用临时梯度值计算获得新的动量项；
利用新的动量项对参数进行调整以获得新的参数值。

基于动量项的参数更新原理参考：
神经网络之反向传播算法（加入动量momentum的误差反向传播算法）
梯度g的计算原理参考：
神经网络之反向传播算法（梯度、误差反向传播算法BP）

2、结合Nesterov动量的误差反向传播算法

将加入动量法的误差反向传播算法应用于神经网络参数优化时的算法步骤如下：

随机初始化神经网络中所有的参数；
设置学习率、动量参数，设置变量v的初始值为0；
输入数据，利用变量v（动量项）对参数进行临时更新；
按照前向传播过程计算网络的输出；
使用损失函数计算误差值；
根据误差值计算出隐含层、输出层每个参数的临时梯度值；
利用临时梯度值计算变量v；
利用变量v对参数进行更新；
重复步骤3到步骤8，当满足停止迭代条件时输出训练后的参数。

误差反向传播算法BP原理参考：
神经网络之反向传播算法（梯度、误差反向传播算法BP）
参数初始化方法参考：
神经网络基础知识之参数初始化

3、算法实现

以数据预测为例，下面介绍加入Nesterov动量的误差反向传播算法的实现过程，将加入Nesterov动量的误差反向传播算法应用于普通的三层神经网络（输入层、隐含层、输出层）的反向传播过程。
选用某省市的表层土壤重金属元素数据集作为实验数据，该数据集总共96组，随机选择其中的24组作为测试数据集，72组作为训练数据集。选取重金属Ti的含量作为待预测的输出特征，选取重金属Co、Cr、Mg、Pb作为模型的输入特征。

3.1 训练过程

#库的导入
import numpy as np
import pandas as pd

#激活函数tanh
def tanh(x):
    return (np.exp(x)-np.exp(-x))/(np.exp(x)+np.exp(-x))
#激活函数偏导数
def de_tanh(x):
    return (1-x**2)

# 梯度加权平均函数，0.9为衰减系数
def accumulation(v,delta):
    v = 0.9 * v + 0.1 * delta
    return  v

# 参数更新量计算函数
def adjust(v):
    change = (-0.01)*v
    return change

#输入数据的导入
df = pd.read_csv("train.csv")
df.columns = ["Co", "Cr", "Mg", "Pb", "Ti"]
Co = df["Co"]
Co = np.array(Co)
Cr = df["Cr"]
Cr = np.array(Cr)
Mg=df["Mg"]
Mg=np.array(Mg)
Pb = df["Pb"]
Pb =np.array(Pb)
Ti = df["Ti"]
Ti = np.array(Ti)
samplein = np.mat([Co,Cr,Mg,Pb])

#数据归一化，将输入数据压缩至0到1之间，便于计算，后续通过反归一化恢复原始值
sampleinminmax = np.array([samplein.min(axis=1).T.tolist()[0],samplein.max(axis=1).T.tolist()[0]]).transpose()
sampleout = np.mat([Ti])
sampleoutminmax = np.array([sampleout.min(axis=1).T.tolist()[0],sampleout.max(axis=1).T.tolist()[0]]).transpose()
sampleinnorm = (2*(np.array(samplein.T)-sampleinminmax.transpose()[0])/(sampleinminmax.transpose()[1]-sampleinminmax.transpose()[0])-1).transpose()
sampleoutnorm = (2*(np.array(sampleout.T)-sampleoutminmax.transpose()[0])/(sampleoutminmax.transpose()[1]-sampleoutminmax.transpose()[0])-1).transpose()
sampleinmax = np.array([sampleinnorm.max(axis=1).T.tolist()]).transpose()
sampleinmin = np.array([sampleinnorm.min(axis=1).T.tolist()]).transpose()
#为归一化后的数据添加噪声
noise = 0.03*np.random.rand(sampleoutnorm.shape[0],sampleoutnorm.shape[1])
sampleoutnorm += noise
sampleinnorm = np.mat(sampleinnorm)

maxepochs = 1000  #训练次数
errorfinal = 0.65*10**(-3)  #停止训练误差阈值
samnum = 72  #输入数据数量
indim = 4  #输入层节点数
outdim = 1  #输出层节点数
hiddenunitnum = 8  #隐含层节点数

#利用归一化后的输入数据初始化参数w1、b1、w2、b2
dvalue = sampleinmax-sampleinmin
valuemid=(sampleinmin+sampleinmax)/2
wmag=0.7*(hiddenunitnum**(1/indim))
rand1=np.random.rand(hiddenunitnum,outdim)
rand2=np.random.randn(hiddenunitnum,indim)
rand1=rand1*wmag
rand2=rand2*wmag
b1=rand1-np.dot(rand2,valuemid)
for i in range(hiddenunitnum):
    for j in range(indim):
        rand2[i][j]=(2*rand2[i][j])/dvalue[j]
w1=rand2
w2 = np.random.uniform(low=-1, high=1, size=[outdim,hiddenunitnum])
b2 = np.random.uniform(low=-1, high=1, size=[outdim,1])

#参数w1、b1、w2、b2均为矩阵形式参与计算，其形状依次为8*4，8*1，1*8，1*1
w1 = np.mat(w1)
b1 = np.mat(b1)
w2 = np.mat(w2)
b2 = np.mat(b2)

#errhistory存储每次训练后的预测值与真实值的误差
errhistory = []

#vw2、vb2、vw1、vb1分别保存参数w1、b1、w2、b2的梯度加权平均值，其形状与w1、b1、w2、b2一一对应
vw2 = np.zeros((1,8))
vb2 = np.zeros((1,1))
vw1 = np.zeros((8,4))
vb1 = np.zeros((8,1))

#deltaw1、deltab1、deltaw2 、deltab2分别保存参数w1、b1、w2、b2的临时更新量
deltaw2 = np.zeros((1,8))
deltab2 = np.zeros((1,1))
deltaw1 = np.zeros((8,4))
deltab1 = np.zeros((8,1))

for i in range(maxepochs):
    #前向传播
    # 利用参数临时更新量对参数w1、b1、w2、b2进行临时更新
    w1 += deltaw1
    b1 += deltab1
    w2 += deltaw2
    b2 += deltab2
    #计算隐含层输出hiddenout，输出层输出networkout
    hiddenout = tanh((np.dot(w1,sampleinnorm).transpose()+b1.transpose())).transpose()
    networkout = np.dot(w2,hiddenout).transpose()+b2.transpose()
    for j in range(samnum):
        networkout[j,:] = tanh(networkout[j,:])
    networkout = networkout.transpose()
    # 计算损失函数
    err = sampleoutnorm - networkout
    loss = np.sum(np.abs(err))/samnum
    sse = np.sum(np.square(err))
    #判断是否满足停止训练条件
    errhistory.append(sse)
    if sse < errorfinal:
        break
    #反向传播
    #利用损失函数计算结果和激活函数偏导数，来计算参数w1、b1、w2、b2的梯度值
    delta2 = np.zeros((outdim,samnum))
    for n in range(samnum):
        delta2[:,n] = (-1) * err[:,n] * de_tanh(networkout[:,n])
    delta1 = np.zeros((hiddenunitnum,samnum))
    for e in range(samnum):
        for f in range(hiddenunitnum):
            delta1[f,e] = w2[:,f] * delta2[:,e] * de_tanh(hiddenout[f,e])
    dw2now = np.dot(delta2,hiddenout.transpose()) #1*8
    db2now = np.dot(delta2,np.ones((samnum,1))) #1*1
    dw1now = np.dot(delta1,sampleinnorm.transpose()) #8*4
    db1now = np.dot(delta1,np.ones((samnum,1))) #8*1

    # 先更新输出层参数
    # w2更新，依次更新w2的梯度加权平均值，参数更新量，w2
    for m in range(hiddenunitnum):
        vw2[:,m] = accumulation(vw2[:,m],dw2now[:,m])
        deltaw2[:,m] = adjust(vw2[:,m])
        w2[:,m] = w2[:,m] + deltaw2[:,m]

    #b2更新，依次更新b2的梯度加权平均值，参数更新量，b2
    vb2 = accumulation(vb2,db2now)
    deltab2 = adjust(vb2)
    b2 = b2 + deltab2

    # 更新隐含层参数
    #w1更新，依次更新w1的梯度加权平均值，参数更新量，w1
    for a in range(hiddenunitnum):
        for b in range(indim):
            vw1[a,b] = accumulation(vw1[a,b],dw1now[a,b])
            deltaw1[a,b] = adjust(vw1[a,b])
            w1[a,b] = w1[a,b] + deltaw1[a,b]

    #b1更新，依次更新b1的梯度加权平均值，参数更新量，b1
    for n in range(hiddenunitnum):
        vb1[n,:] = accumulation(vb1[n,:],db1now[n,:])
        deltab1[n,:] = adjust(vb1[n,:])
        b1[n,:] = b1[n,:] + deltab1[n,:]


    print("the generation is:",i,",the loss is:",loss)

#达到最大训练次数，保存此时的参数w1、b1、w2、b2
np.save("w1.npy",w1)
np.save("b1.npy",b1)
np.save("w2.npy",w2)
np.save("b2.npy",b2)