【机器学习】梯度上升算法实验

2020-10-24

逻辑回归迭代过程权重变化绘图

import numpy as np
from matplotlib import pyplot as plt

def sigmoid(inX):
    return 1.0 / (1 + np.exp(-inX))

def loadDataSet():
dataMat = []; labelMat = []
fr = open('testSet.txt')
for line in fr.readlines():
    lineArr = line.strip().split()
    dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
    labelMat.append(int(lineArr[2]))
return dataMat, labelMat


def plotBestFit(weights):
    dataMat, labelMat = loadDataSet()
    dataArr = np.array(dataMat)
    n = np.shape(dataArr)[0]
    xcord1 = []; ycord1 = []
    xcord2 = []; ycord2 = []
    for i in range(n):
        if int(labelMat[i]) == 1:
            xcord1.append(dataArr[i, 1]); ycord1.append(dataArr[i, 2])
        else:
            xcord2.append(dataArr[i, 1]); ycord2.append(dataArr[i, 2])
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
    ax.scatter(xcord2, ycord2, s=30, c='green')
    x = np.arange(-3.0, 3.0, 0.1)
    y = (-weights[0]-weights[1]*x)/weights[2]  #得到回归曲线
    ax.plot(x, y)
    plt.xlabel('X1'); plt.ylabel('X2')
    plt.show()

梯度下降算法

def gradAscent(dataMatIn, classLabels,maxCycles=500):
    dataMatrix = np.mat(dataMatIn)             #convert to NumPy matrix
    labelMat = np.mat(classLabels).transpose() #convert to NumPy matrix
    m, n = np.shape(dataMatrix)
    alpha = 0.001
    weightlist = []
    weights = np.ones((n, 1))
    for k in range(maxCycles):              #heavy on matrix operations
        weightlist.append(weights)
        h = sigmoid(dataMatrix*weights)     #matrix mult
        error = (labelMat - h)              #vector subtraction
        weights = weights + alpha * dataMatrix.transpose()* error #matrix mult
    return np.array(weightlist)


def stocGradAscent0(dataMatrix, classLabels,batch = 1500):
    dataMatrix = np.array(dataMatrix)#转换为numpy矩阵
    m, n = np.shape(dataMatrix)#得到矩阵维度
    alpha = 0.01#学习率
    weights = np.ones(n)   #initialize to all ones  1x3
    weightlist = []#权重列表
    for _ in range(batch):#重复迭代
        for i in range(m):
            weightlist.append(weights)
            h = sigmoid(np.sum(dataMatrix[i]*weights))#得到预测值
            error = classLabels[i] - h# 得到与目标的误差值 1x1
            weights = weights + alpha * error * dataMatrix[i]#反向传播，更新权重值
    return np.array(weightlist)#转换为numpy矩阵返回

def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    dataMatrix = np.array(dataMatrix)
    m, n = np.shape(dataMatrix)
    weights = np.ones(n)   #initialize to all ones
    weightlist = []
    for j in range(numIter):
        dataIndex = list(range(m))
        for i in range(m):
            weightlist.append(weights)
            alpha = 4/(1.0+j+i)+0.0001    #apha decreases with iteration, does not
            randIndex = int(np.random.uniform(0, len(dataIndex)))#go to 0 because of the constant
            h = sigmoid(sum(dataMatrix[randIndex]*weights))
            error = classLabels[randIndex] - h
            weights = weights + alpha * error * dataMatrix[randIndex]
            del(dataIndex[randIndex])#删除使用过的数据
    return np.array(weightlist)

逻辑回归随机梯度上升算法权重绘制

绘图算法原理:

在梯度下降算法的开始定义一个列表weightlist = []用于存放每次迭代过程中的==权重==变化，每次迭代使用weightlist.append(weights)将该次迭代的权重值weights存入列表，这样便得到了权重历史迭代数据
绘图;

因为有三个权重，所以使用plt.subplots(nrows=3, ncols=1, figsize=(20, 15))将图分为3x1
X轴为迭代次数，从weightlist长度可以得到X = range(len(weightlist))
Y轴即为迭代的权重数据，W[:, n]得到第n列数据，W.shape = (100,3)

梯度上升算法原理:

梯度上升法的思想：要找到某函数的最大值 ,最好的方法是沿着该函数的梯度方向探寻。

梯度上升算法的迭代公式：$w := w+\alpha*dw (f(w)) $(dw是关于系数w的梯度，$\alpha$是学习率)

每个回归系数初始化为1
重复R次：
       计算整个数据集的梯度
       使用alpha X gradient更新回归系数
返回回归系数

运行结果:

可见两种随机梯度上升算法的权重值W[0],W[1]最后都趋于稳定值。

不同的是没有随机选取样本数据而是顺序遍历每个数据点的stocGradAscent0方法迭代200(每次遍历全部100给数据点)次后任然没有达到足够的稳定，而使用随机算法选取每次迭代数据点的stocGradAscent1方法在约25次的迭代后就开始收敛了。可见改进后的随机梯度下降算法效果较好。

stocGradAscent1效果较好的原因：

第一点增加了alpha动态减少的机制，这样做的原因是为了保证在多次迭代之后新数据仍然具有一定的影响。

第二点是通过随机选取样本来更新回归系数，这种方法减少周期性的波动。每次随机从列表中选出一个值,然后从列表中删掉该值，重新迭代。

#W权重参数列表W.shape = (m,3),title标题
def plotWeightlist2(W,title):
    font ={'family' : 'simsun',
            'weight' : 'bold',
            'size'   : 20,
            }#字体
    m,n = W.shape
    X = range(m)#生成X轴坐标点
    # plt.figure(1)#标号
    plt.subplots(nrows=3, ncols=1, figsize=(20, 15))#分成3x1
    plt.subplot(311)#分成3x1，占用第一个，即第一行第一列的子图
    plt.title(title,font)
    plt.plot(X, W[:, 0])
    plt.ylabel(u'X0', font)
    plt.subplot(312)
    plt.plot(X, W[:, 1], 'r--')
    plt.ylabel(u'X1', font)
    plt.subplot(313)
    plt.plot(X, W[:, 2], )
    plt.xlabel(r'迭代次数', font)
    plt.ylabel(u'X2', font)
    plt.show()

测试

dataset,labelset = loadDataSet()
print("dataset:"+str(len(dataset)))
weightlist = stocGradAscent0(dataset,labelset,200)
plotWeightlist2(weightlist,"随机梯度上升算法")
plotBestFit(weightlist[len(weightlist)-1])#取出最后一次的的权重最为最好的权重绘制回归曲线

weightlist = stocGradAscent1(dataset,labelset,200)
print(weightlist.shape)
plotWeightlist2(weightlist,"改进后的随机梯度上升算法")
plotBestFit(weightlist[len(weightlist)-1])#取出最后一次的的权重最为最好的权重绘制回归曲线

简单神经网络迭代过程

绘图算法原理:

与之前逻辑回归一样，在迭代开始定义一个列表losslist = []用于存放每次迭代过程中的==loss==变化，每次迭代使用losslist.append(loss)将该次迭代的Loss值loss存入列表，这样便得到了Loss值历史迭代数据
绘图与之前一样，不多赘述

梯度下降算法原理:

与梯度上升类似，只不过是计算使得Loss值最小的权重参数，例如$\frac{\partial Loss}{\partial W_1} = 1$说明W1对Loss的贡献是$+1$,那么为了使得Loss值减小：
$W_1 := W_1 - \frac{\partial Loss}{\partial W_1}$,当然这样学习率就默认为1了，可能会离理想的最小Loss值点较远，所以需要设置学习率$\alpha$,做最终就是$W_1 := W_1 - \alpha * \frac{\partial Loss}{\partial W_1}$。实际过程中可能会加入例如ReLU,Sigmoid等激活函数，包含读个隐含层等，这样使用链式求导法则更新各级权重(反向传播)即可。

运行结果:
可见随着迭代次数的增加，Loss值不断减小，符合预期。

import numpy as np
from matplotlib import pyplot as plt

def sigmoid(x):
  # Sigmoid activation function: f(x) = 1 / (1 + e^(-x))
  return 1 / (1 + np.exp(-x))

def deriv_sigmoid(x):
  # Derivative of sigmoid: f'(x) = f(x) * (1 - f(x))
  fx = sigmoid(x)
  return fx * (1 - fx)

def mse_loss(y_true, y_pred):
  # y_true and y_pred are numpy arrays of the same length.
  return ((y_true - y_pred) ** 2).mean()

class OurNeuralNetwork:
  '''
  A neural network with:
    - 2 inputs
    - a hidden layer with 2 neurons (h1, h2)
    - an output layer with 1 neuron (o1)

  *** DISCLAIMER ***:
  The code below is intended to be simple and educational, NOT optimal.
  Real neural net code looks nothing like this. DO NOT use this code.
  Instead, read/run it to understand how this specific network works.
  '''
  def __init__(self):
    # Weights
    self.w1 = np.random.normal()
    self.w2 = np.random.normal()
    self.w3 = np.random.normal()
    self.w4 = np.random.normal()
    self.w5 = np.random.normal()
    self.w6 = np.random.normal()

    # Biases
    self.b1 = np.random.normal()
    self.b2 = np.random.normal()
    self.b3 = np.random.normal()

  def feedforward(self, x):
    # x is a numpy array with 2 elements.
    h1 = sigmoid(self.w1 * x[0] + self.w2 * x[1] + self.b1)
    h2 = sigmoid(self.w3 * x[0] + self.w4 * x[1] + self.b2)
    o1 = sigmoid(self.w5 * h1 + self.w6 * h2 + self.b3)
    return o1

  def train(self, data, all_y_trues):
    '''
    - data is a (n x 2) numpy array, n = # of samples in the dataset.
    - all_y_trues is a numpy array with n elements.
      Elements in all_y_trues correspond to those in data.
    '''
    learn_rate = 0.1
    epochs = 1000 # number of times to loop through the entire dataset

    losslist = []
    for epoch in range(epochs):
      for x, y_true in zip(data, all_y_trues):
        # --- Do a feedforward (we'll need these values later)
        sum_h1 = self.w1 * x[0] + self.w2 * x[1] + self.b1
        h1 = sigmoid(sum_h1)

        sum_h2 = self.w3 * x[0] + self.w4 * x[1] + self.b2
        h2 = sigmoid(sum_h2)

        sum_o1 = self.w5 * h1 + self.w6 * h2 + self.b3
        o1 = sigmoid(sum_o1)
        y_pred = o1

        # --- Calculate partial derivatives.
        # --- Naming: d_L_d_w1 represents "partial L / partial w1"
        d_L_d_ypred = -2 * (y_true - y_pred)

        # Neuron o1
        d_ypred_d_w5 = h1 * deriv_sigmoid(sum_o1)
        d_ypred_d_w6 = h2 * deriv_sigmoid(sum_o1)
        d_ypred_d_b3 = deriv_sigmoid(sum_o1)

        d_ypred_d_h1 = self.w5 * deriv_sigmoid(sum_o1)
        d_ypred_d_h2 = self.w6 * deriv_sigmoid(sum_o1)

        # Neuron h1
        d_h1_d_w1 = x[0] * deriv_sigmoid(sum_h1)
        d_h1_d_w2 = x[1] * deriv_sigmoid(sum_h1)
        d_h1_d_b1 = deriv_sigmoid(sum_h1)

        # Neuron h2
        d_h2_d_w3 = x[0] * deriv_sigmoid(sum_h2)
        d_h2_d_w4 = x[1] * deriv_sigmoid(sum_h2)
        d_h2_d_b2 = deriv_sigmoid(sum_h2)

        # --- Update weights and biases
        # Neuron h1
        self.w1 -= learn_rate * d_L_d_ypred * d_ypred_d_h1 * d_h1_d_w1
        self.w2 -= learn_rate * d_L_d_ypred * d_ypred_d_h1 * d_h1_d_w2
        self.b1 -= learn_rate * d_L_d_ypred * d_ypred_d_h1 * d_h1_d_b1

        # Neuron h2
        self.w3 -= learn_rate * d_L_d_ypred * d_ypred_d_h2 * d_h2_d_w3
        self.w4 -= learn_rate * d_L_d_ypred * d_ypred_d_h2 * d_h2_d_w4
        self.b2 -= learn_rate * d_L_d_ypred * d_ypred_d_h2 * d_h2_d_b2

        # Neuron o1
        self.w5 -= learn_rate * d_L_d_ypred * d_ypred_d_w5
        self.w6 -= learn_rate * d_L_d_ypred * d_ypred_d_w6
        self.b3 -= learn_rate * d_L_d_ypred * d_ypred_d_b3

      # --- Calculate total loss at the end of each epoch
        y_preds = np.apply_along_axis(self.feedforward, 1, data)
        loss = mse_loss(all_y_trues, y_preds)
        losslist.append(loss)
      if epoch % 100 == 0:
        print("Epoch %d loss: %.3f" % (epoch, loss))
    return np.array(losslist)
# Define dataset
data = np.array([
  [-2, -1],  # Alice
  [25, 6],   # Bob
  [17, 4],   # Charlie
  [-15, -6], # Diana
])
all_y_trues = np.array([
  1, # Alice
  0, # Bob
  0, # Charlie
  1, # Diana
])


def plotLosslist(Loss, title):
  font = {'family': 'simsun',
          'weight': 'bold',
          'size': 20,
          }
  m = len(Loss)
  X = range(m)
  # plt.figure(1)
  plt.subplots(nrows=1, ncols=1, figsize=(10, 8))
  plt.subplot(111)
  plt.title(title, font)
  plt.plot(X, Loss)
  plt.xlabel(r'迭代次数', font)
  plt.ylabel(u'Loss', font)
  plt.show()

def drawweightlist(X, title):
  X = X.reshape(X.shape[0])
  plt.title(title)
  plt.plot(range(X.shape[0]), X)
  plt.show()

# Train our neural network!
network = OurNeuralNetwork()
losslist = network.train(data, all_y_trues)

plotLosslist(losslist,"Loss迭代图")
# Make some predictions
emily = np.array([-7, -3]) # 128 pounds, 63 inches
frank = np.array([20, 2])  # 155 pounds, 68 inches
print("Emily: %.3f" % network.feedforward(emily)) # 0.951 - F
print("Frank: %.3f" % network.feedforward(frank)) # 0.039 - M

关于《机器学习实践-Peter Harrington》

我认为书上随机梯度上升算法的迭代图的第一幅图是错误的！！

理由如下：
书上给出的回归曲线y = (-weights[0]-weights[1]*x)/weights[2]，显然对于回归曲线$\hat{y} = \hat{a} x+b$,$\hat{a}$应该是相近的，至少正负号应该是一样的。按照书上给出的代码 $\hat{a}$ = -weights[1]//weights[2]，观察发现$\hat{a} < 0$。这就非常奇怪的，所以我认为是书的错误。

| 阅次 | 字数:2,603