【python】随机产生20个点用梯度下降法线性回归拟合---BGD+SGD+MBGD

2021/04/07 14:57:42

目录

  • 什么是梯度下降法
  • 怎么用梯度下降法进行拟合(以BGD为例)
  • 其他改进形式梯度下降法(SGD+MBGD)

什么是梯度下降法

梯度下降算法原理讲解——机器学习
原理网上有很多,这个博客比较详细友好

怎么用梯度下降法进行拟合(以BGD为例)

一道作业题:
随机产生20个点,用线性回归拟合,并画出迭代次数与总损失值的关系曲线图和拟合结果图。

  1. 怎么拟合一道直线呢?先把直线方程设出来,h为预测函数
    在这里插入图片描述
    现在需要求解最佳的 θ 0 \theta_0 θ0 θ 1 \theta_1 θ1得到最佳的直线

  2. 什么是最佳直线,就是用预测函数 h h h算出的预测值与实际值的差之和最小(只考虑y值), J J J代价函数
    在这里插入图片描述
    m是数据集中数据点的个数,也就是样本数;
    ½是一个常量,这样是为了在求梯度的时候,二次方乘下来的2就和这里的½抵消了,自然就没有多余的常数系数,方便后续的计算,同时对结果不会有影响。

  3. 那么现在就是求上述代价函数 J J J(二元函数)的最小值,怎么求呢?用梯度下降法求,求出最佳的 θ 0 \theta_0 θ0 θ 1 \theta_1 θ1使得代价函数 J J J取得最小值。按照如下的迭代公式求解(注意 θ \theta θ的上标表示先后两次迭代)
    在这里插入图片描述
    α \alpha α为学习率or步长,定义每一步的长度,不能太短也不能太长

  4. 为了简化计算,我们将 J J J ∇ J \nabla J J矩阵形式表示,如下
    在这里插入图片描述在这里插入图片描述
    如何随机产生点

# 20个点
m = 20
# 构造矩阵
x0 = np.ones((m,1))
x1 = np.arange(1, m+1).reshape(m, 1)
x = np.hstack((x0, x1)) # 得到x的值
y = np.random.randint(0, 30, size=(m, 1))
y.sort(axis=0) # 得到随机产生y的值,排序是为了使其类似直线分布

效果如图:
在这里插入图片描述
代价函数

# 求损失值函数
def loss_function(p, x, y):
    temp = np.dot(x, p) - y
    result = 1 / (2*m) * np.dot(np.transpose(temp), temp)
    return result

计算梯度函数

# 求梯度函数
def grad_function(p, x, y):
    temp = np.dot(x, p) - y
    result = 1 / m * np.dot(np.transpose(x),temp)
    return result

梯度迭代过程

# 存储迭代次数
count = 0 
# 设置p的初始值
p = np.array([1, 1]).reshape(2, 1)
# 计算初始梯度
grad = grad_function(p, x, y)
# 计算初始损失值
loss = loss_function(p, x, y)
# 存储损失值
lloss = []
loss = loss_function(p,x,y).flatten().tolist() # 先扁平化得到一维矩阵再转换为列表
lloss.extend(loss) # 这里用extend不用append
count += 1
# 设定精度
e = 1e-5
# 开始循环迭代
while np.all(np.absolute(grad) > e):
    p = p - algha * grad # 迭代公式
    # 更新grad和loss
    grad = grad_function(p, x, y)
    loss = loss_function(p, x, y)
    loss = loss_function(p,x,y).flatten().tolist()
    lloss.extend(loss)
    count += 1

完整代码如下:

import numpy as np 
import matplotlib.pyplot as plt
from numpy.core.numeric import ones 
# 显示中文
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 20个点
m = 20
# 构造矩阵
x0 = np.ones((m,1))
x1 = np.arange(1, m+1).reshape(m, 1)
x = np.hstack((x0, x1)) # 得到x的值
y = np.random.randint(0, 30, size=(m, 1))
y.sort(axis=0) # 得到随机产生y的值,排序是为了使其类似直线分布
# 设定学习率
algha = 0.01
# 求损失值函数
def loss_function(p, x, y):
    temp = np.dot(x, p) - y
    result = 1 / (2*m) * np.dot(np.transpose(temp), temp)
    return result
# 求梯度函数
def grad_function(p, x, y):
    temp = np.dot(x, p) - y
    result = 1 / m * np.dot(np.transpose(x),temp)
    return result

# 存储迭代次数
count = 0 
# 设置p的初始值
p = np.array([1, 1]).reshape(2, 1)
# 计算初始梯度
grad = grad_function(p, x, y)
# 计算初始损失值
loss = loss_function(p, x, y)
# 存储损失值
lloss = []
loss = loss_function(p,x,y).flatten().tolist() # 先扁平化得到一维矩阵再转换为列表
lloss.extend(loss) # 这里用extend不用append
count += 1
# 设定精度
e = 1e-5
# 开始循环迭代
while np.all(np.absolute(grad) > e):
    p = p - algha * grad # 迭代公式
    # 更新grad和loss
    grad = grad_function(p, x, y)
    loss = loss_function(p, x, y)
    loss = loss_function(p,x,y).flatten().tolist()
    lloss.extend(loss)
    count += 1
# 绘图
fig,ax = plt.subplots(1,2) # 一行两列的子图
fig.suptitle('基于梯度下降算法的线性回归拟合')
fig.dpi = 150
# 得到最佳的p
[p0, p1] = p
print('p:', p0,p1)
# 绘制拟合图像
ax[0].set_title('拟合图像')
ax[0].scatter(x1,y)
ax[0].plot(x1, p0 + p1*x1, color='r')
ax[0].grid(True) # 绘制网格线
# 绘制损失值图像
ax[1].set_title('损失值图像')
li = list(range(count))
print('迭代次数:', count)
ax[1].plot(li,lloss)
ax[1].grid(True)
plt.show()

效果如下:
在这里插入图片描述
上述代码中numpy库函数调用的一些解释
Python numpy.transpose 详解
08_Numpy初始化生成相同元素值的ndarray数组
np.dot()函数的用法详解

其他改进形式梯度下降法(SGD+MBGD)

用关键段代码进行一个比较
BGD非矩阵形式
每次使用全量的样本数据进行计算,得到m次梯度的和,用这个和乘以步长进行迭代
count+1是计算一次m个梯度的和,用和乘以步长更新一次 θ \theta θ
![在这里插入图片描述](https://img-blog.csdnimg.cn/20210404205201136.png

代码关键段:(完整附在后面)

while count < loop_max:
    count += 1
    # 在标准梯度下降中,权值更新的每一步对多个样例求和,需要更多的计算
    sum_m = np.zeros(2)
    for i in range(m):
        diff = (np.dot(p, x[i]) - y[i]) * x[i] # 一维向量进行内积,就是对应元素相乘再相加,二维及以上就是矩阵乘法运算
        sum_m = sum_m + diff     # 当alpha取值过大时,sum_m会在迭代过程中会溢出

    p = p - alpha * sum_m       # 注意步长alpha的取值,过大会导致振荡
    #w = w - 0.005 * sum_m # alpha取0.005时产生振荡,需要将alpha调小

    # 判断是否已收敛,前后两次w的值的差小于某个阈值
    if np.linalg.norm(p - error) < epsilon:
        break
    else:
        error = p

SGD
每次从训练集中随机选择一个样本计算得到梯度,用该梯度乘以步长进行迭代
count+1是计算一个数据的梯度,用梯度乘以步长更新一次 θ \theta θ,这个过程进行m次,即总共更新m次 θ \theta θ
在这里插入图片描述

while count < loop_max:
    count += 1
    # 遍历训练数据集,不断更新p
    for i in range(m):
        # 训练集代入,计算误差值
        diff = (np.dot(p, x[i]) - y[i]) * x[i]  
        # 采用随机梯度下降算法,更新一次权值只使用一个训练数据
        p = p - alpha * diff
    # 终止条件:前后两次计算出的权向量的绝对误差充分小
    if np.linalg.norm(p - error) < epsilon:     
        break
    else:
        error = p

MBGD
综合了BGD与SGD,在每次更新速度与更新次数中间取得一个平衡,其每次更新从训练集中随机选择b,b<m个样本计算梯度和,梯度和乘以步长进行迭代。
count+1是计算b个数据的梯度和,用1/b乘以梯度和乘以步长更新 θ \theta θ。这个过程进行m/b次,即更新m/b次 θ \theta θ
在这里插入图片描述

while count < loop_max:
    count += 1
    for i in range(1,m,minibatch_size):
        sum_m = np.zeros(2)
        for k in range(i-1,i+minibatch_size-1,1):
            diff = (np.dot(p, x[k]) - y[k]) *x[k]
            sum_m = sum_m + diff  #当alpha取值过大时,sum_m会在迭代过程中会溢出
        p = p- alpha * (1.0/minibatch_size) * sum_m #注意步长alpha的取值,过大会导致振荡

    # 判断是否已收敛
    if np.linalg.norm(p- error) < epsilon:
        break
    else:
        error = p

完整BGD代码:

# 非矩阵形式实现
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

# 构造训练数据
x1 = np.arange(0., 10., 0.2)
m = len(x1)                                      # 训练数据点数目
x0 = np.full(m, 1.0)
x = np.vstack([x0, x1]).T               # 将偏置b作为权向量的第一个分量
y = 2 * x1 + 5 + np.random.randn(m)

# 两种终止条件
loop_max = 100000   # 最大迭代次数(防止死循环)
epsilon = 1e-5

# 初始化权值
np.random.seed(0)
p = np.random.randn(2)

alpha = 0.001      # 步长(注意取值过大会导致振荡,过小收敛速度变慢)
diff = 0.
error = np.zeros(2)
count = 0          # 循环次数

while count < loop_max:
    count += 1
    # 在标准梯度下降中,权值更新的每一步对多个样例求和,需要更多的计算
    sum_m = np.zeros(2)
    for i in range(m):
        diff = (np.dot(p, x[i]) - y[i]) * x[i] # 一维向量进行内积,就是对应元素相乘再相加,二维及以上就是矩阵乘法运算
        sum_m = sum_m + diff     # 当alpha取值过大时,sum_m会在迭代过程中会溢出

    p = p - alpha * sum_m       # 注意步长alpha的取值,过大会导致振荡
    #w = w - 0.005 * sum_m # alpha取0.005时产生振荡,需要将alpha调小

    # 判断是否已收敛,前后两次w的值的差小于某个阈值
    if np.linalg.norm(p - error) < epsilon:
        break
    else:
        error = p
print ('loop count = %d' % count,  '\tw:[%f, %f]' % (p[0], p[1]))

plt.scatter(x1, y)
plt.plot(x1, p[1] * x1 + p[0], 'r')
plt.show()

完整SGD代码:

import numpy as np
import matplotlib.pyplot as plt

# 构造训练数据(另一种构造方法)
# x1 = np.arange(0.,10.,0.2)
# m = len(x1) # 训练数据点数目
# x0 = np.full(m, 1.0)
# x = np.vstack([x0, x1]).T # 将偏置b作为权向量的第一个分量
# y = 2 *x + 5 +np.random.randn(m) # 基于函数y=2x+5

m = 20
x0 = np.full(m, 1.0)
x1 = np.arange(1, m+1)  
x = np.vstack([x0, x1]).T  
y = np.random.randint(0, 30,size=20)
y.sort(axis=0)

# 两种终止条件
loop_max = 10000   # 最大迭代次数(防止死循环)
epsilon = 1e-5

# 设置p的初始值
np.random.seed(0)
p = np.random.randn(2)
# 设定学习率(注意取值过大会导致振荡,过小收敛速度变慢)
alpha = 0.001   
# 精度 
error = np.zeros(2)
# 循环次数
count = 0         
# 开始迭代
while count < loop_max:
    count += 1
    # 遍历训练数据集,不断更新p
    for i in range(m):
        # 训练集代入,计算误差值
        diff = (np.dot(p, x[i]) - y[i]) * x[i]  
        # 采用随机梯度下降算法,更新一次权值只使用一个训练数据
        p = p - alpha * diff
    # 终止条件:前后两次计算出的权向量的绝对误差充分小
    if np.linalg.norm(p - error) < epsilon:     
        break
    else:
        error = p
[p0, p1] = p
print('p:', p0,p1)
print('迭代次数:', count)
# 绘图
plt.scatter(x1, y)
plt.plot(x1,  p0 + p1 * x1, 'r')
plt.show()

完整MBGD代码

import numpy as np
import matplotlib.pyplot as plt

# 构造训练数据
# x1 = np.arange(0.,10.,0.2)
# m = len(x1) # 训练数据点数目
# x0 = np.full(m, 1.0)
# x = np.vstack([x0, x1]).T # 将偏置b作为权向量的第一个分量
# y = 2 * x1 + 5 +np.random.randn(m)

m = 20
x0 = np.full(m, 1)
x1 = np.arange(1, m+1)  
x = np.vstack([x0, x1]).T  
y = np.random.randint(0, 30,size=20)
y.sort(axis=0)

# 两种终止条件
loop_max = 10000  #最大迭代次数(防止死循环)
epsilon = 1e-5

# 设置p的初始值
np.random.seed(0)
p = np.random.randn(2)
#步长(注意取值过大会导致振荡即不收敛,过小收敛速度变慢)
alpha = 0.001  
diff = 0.
error = np.zeros(2)
# 循环次数
count = 0  
# 每次更新的样本数
minibatch_size = 5  
while count < loop_max:
    count += 1
    # 在minibatch梯度下降中,权值更新的每一步对多个样例求和,需要更多的计算
    for i in range(1,m,minibatch_size):
        sum_m = np.zeros(2)
        for k in range(i-1,i+minibatch_size-1,1):
            diff = (np.dot(p, x[k]) - y[k]) *x[k]
            sum_m = sum_m + diff  #当alpha取值过大时,sum_m会在迭代过程中会溢出
        p = p- alpha * (1.0/minibatch_size) * sum_m #注意步长alpha的取值,过大会导致振荡

    # 判断是否已收敛
    if np.linalg.norm(p- error) < epsilon:
        break
    else:
        error = p

[p0, p1] = p
print('p:', p0,p1)
print('迭代次数:', count)
# 绘制
plt.scatter(x1, y)
plt.plot(x1, p[1]* x1 + p[0],'r')
plt.show()

参考梯度下降算法及python实现(学习笔记)