梯度下降法(Gradient descent)_def gradient_descent(w, b, x, y, eta):-程序员宅基地

1. 梯度下降法简介

1-1

 

以下是定义了一个损失函数以后,参数theta对应的损失函数J的值对应的示例图,我们需要找到使得损失函数值J取得最小值对应的theta(这里是二维平面,也就是我们的参数只有一个)

在直线方程中,导数代表斜率
在曲线方程中,导数代表切线斜率
导数代表theta单位变化时,J相应的变化

 

1-2

 

1-3

 

η太小,会减慢收敛学习速度

 

1-4


η太大,甚至导致不收敛

1-5

其他注意事项

  • 并不是所有函数都有唯一的极值点

     

    1-6

     

    解决方案:

    • 多次运行,随机化初始点
    • 梯度下降法的初始点也是一个超参数

1-7


2 梯度下降法模拟

2.1 实现

import numpy as np
import matplotlib.pyplot as plt
# 简单模拟一个损失函数
plot_x = np.linspace(-1,6,141)
plot_y = (plot_x-2.5)**2-1
# 绘制我们模拟的损失函数
plt.plot(plot_x,plot_y)

2.1-1

def dJ(theta):
    """损失函数的导数"""
    return 2*(theta-2.5)

def J(theta):
    """损失函数"""
    try:
        return (theta-2.5)**2-1
    except:
        return float('inf')

def gradient_descent(initial_theta,eta,n_iters = 1e4,epsilon=1e-8):
    """
    梯度下降法封装
    initial_theta:初始化的theta值
    eta:学习率η
    n_iters: 最大循环次数
    epsilon: 精度
    """
    theta = initial_theta
    # theta_history 保存theta的变化值
    theta_history.append(initial_theta)
    i_iters = 0
    
    while i_iters<n_iters:
        """
        如果theta两次变化之间的损失函数值的变化小于我们定义的精度
        则可以说明我们已经找到了最低的损失函数值和对应的theta
        
        如果循环次数超过了我们设置的循环次数,
        则说明可能由于η设置的过大导致无止境的循环
        """
        gradient = dJ(theta)
        last_theta = theta
        theta = theta - eta * gradient
        theta_history.append(theta)
        
        if (abs(J(theta)-J(last_theta)) < epsilon):
            break
        i_iters += 1

def plot_theta_history():
    plt.plot(plot_x,J(plot_x))
    plt.plot(np.array(theta_history),J(np.array(theta_history)),color='r',marker='+')
    print("the size of theta_history is %d"%len(theta_history))

2.2 使用不同η学习率测试并观察我们的梯度下降法的结果

eta = 0.1
theta_history = []
gradient_descent(0.,eta)
plot_theta_history()

2.2-1

eta = 0.01
theta_history = []
gradient_descent(0.,eta)
plot_theta_history()

2.2-2

eta = 0.001
theta_history = []
gradient_descent(0.,eta)
plot_theta_history()

2.2-3

eta = 0.8
theta_history = []
gradient_descent(0.,eta)
plot_theta_history()

2.2-4

 

可以发现,只要η不超过一个限度,我们编写的函数都可以在有限次数之后找到最优解,并且η越小,学习的次数越多
下面来看一下,如果eta取较大值1.1,会出现什么情况

eta = 1.1
theta_history = []
gradient_descent(0.,eta)
# 数据量太大会报错
# plot_theta_history()
print(len(theta_history))
# 输出10001
theta_history[-1]
# 输出 nan(not a number)

可以看出当我们的eta取1.1,函数会循环直至终止,这是由于,我们的η设置过大,导致每次循环过后,损失函数j的值都向大的方向变化

eta = 1.1
theta_history = []
gradient_descent(0.,eta,n_iters = 10)
plot_theta_history()

2.2-5


3 多元线性回归中的梯度下降法

3-1

 

一个三维空间中的梯度下降法(x,y为系数,z为损失函数)

 

3-2

 

推导过程

 

3-3

 

3-4


上面推导出的式子的大小是和样本数有关的,m越大,结果越大,这是不合理的,我们希望和m无关

3-5


4 线性回归中的梯度下降法的实现

4.1 比较笨的方法实现

def fit_gd(self, X_train, y_train, eta=0.01, n_iters = 1e4):
        """根据训练数据集X_train,y_train, 使用梯度下降法训练Linear Regression 模型"""
        assert X_train.shape[0] == y_train.shape[0], \
            "the size of X_train must be equal to the size of y_train"

        def J(theta, X_b, y):
            try:
                return np.sum((y - X_b.dot(theta))**2) / len(X_b)
            except:
                return float('inf')

        def dJ(theta, X_b, y):
            res = np.empty(len(theta))
            res[0] = np.sum(X_b.dot(theta) - y)

            for i in range(1, len(theta)):
                res[i] = np.sum((X_b.dot(theta) - y).dot(X_b[:, i]))

            return res * 2 / len(X_b)

        def gradient_descent(X_b, y, initial_theta, eta, n_iters=n_iters, epsilon=1e-8):
            """
            梯度下降法封装
            X_b: X特征矩阵
            y: 结果向量
            initial_theta:初始化的theta值
            eta:学习率η
            n_iters: 最大循环次数
            epsilon: 精度
            """
            theta = initial_theta
            i_iters = 0

            while i_iters < n_iters:
                """
                如果theta两次变化之间的损失函数值的变化小于我们定义的精度
                则可以说明我们已经找到了最低的损失函数值和对应的theta
                
                如果循环次数超过了我们设置的循环次数,
                则说明可能由于η设置的过大导致无止境的循环
                """
                gradient = dJ(theta, X_b, y)
                last_theta = theta
                theta = theta - eta * gradient

                if abs(J(theta, X_b, y) - J(last_theta, X_b, y)) < epsilon:
                    break

                i_iters += 1

            return theta

        X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
        initial_theta = np.zeros(X_b.shape[1])
        self._theta = gradient_descent(X_b, y_train, initial_theta, eta)

        self.interception_ = self._theta[0]
        self.coef_ = self._theta[1:]

        return self

4.2 测试我们的算法

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(666)
x = 2*np.random.random(size = 100)
# 定义截距为4 斜率为3
y = x * 3. + 4. + np.random.normal(size=100)
X = x.reshape(-1,1)
plt.scatter(x,y)

4.2-1

from machine_learning.LinearRegression import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit_gd(X,y)

lin_reg.interception_
# 4.021457858204859
lin_reg.coef_
# array([3.00706277])

4.3 向量化

4.3-1

 

4.3-2

 

4.3-3

修改之前的求导函数

def dJ(theta, X_b, y):
            # res = np.empty(len(theta))
            # res[0] = np.sum(X_b.dot(theta) - y)
            #
            # for i in range(1, len(theta)):
            #     res[i] = np.sum((X_b.dot(theta) - y).dot(X_b[:, i]))
            #
            # return res * 2 / len(X_b)
            return X_b.T.dot(X_b.dot(theta) - y) * 2. / len(X_b)

使用真实的数据测试

 

4.3-4

 

使用真实的数据,调整eta和iters,要么由于eta太小导致无法得出真实的结果,要么由于eta太大导致训练时间加长,这是由于数据的规模在不同的特征上不同,所以我们需要对数据进行归一化

4.4 数据归一化

4.4-1

4.4-2

4.4-3

如果样本数非常多,那么即使使用梯度下降法也会导致速度比较慢,因为在梯度下降法中,每一个样本都要参与运算。这时候需要采用随机梯度下降法,我们将在下一小节进行介绍


5 随机梯度下降法

5.1 随机梯度下降法介绍

5.1

 

批量梯度下降法带来的一个问题是η的值需要设置的比较小,在样本数比较多的时候导致不是速度特别慢,这时候观察随机梯度下降法损失函数的求导公式,可以发现,我们对每一个Xb都做了求和操作,又在最外面除以了m,那么可以考虑将求和和除以m的两个运算约掉,采用每次使用一个随机的Xb

 

5.2

 

5.3

 

由于我们使用的事随机梯度下降法,所以导致我们的最终结果不会像批量梯度下降法一样准确的朝着一个方向运算,而是曲线行下降,这时候我们就希望,越到下面,η值相应减小,事运算次数变多,从而精确计算结果

 

5-4


这里使用了模拟退火的思想

5-5

5.2 随机梯度下降法实现

def dJ_sgd(theta,X_b_i,y_i):
    return X_b_i.T.dot(X_b_i.dot(theta) - y_i) * 2

def sgd(X_b,y,initial_theta,n_iters):
    
    t0 = 5
    t1 = 50
    
    def learning_rate(t):
        return t0 / (t + t1)
    
    theta = initial_theta
    for cur_iter in range(n_iters):
        rand_i = np.random.randint(len(X_b))
        gradient = dJ_sgd(theta,X_b[rand_i],y[rand_i])
        theta = theta - learning_rate(cur_iter) * gradient
    return theta
        %%time
        eta = 0.01
        X_b = np.hstack([np.ones((len(X), 1)), X])
        initial_theta = np.zeros(X_b.shape[1])
        # 随机的检查了3分之一个样本总量的样本
        _theta = sgd(X_b, y, initial_theta, n_iters=len(X_b)//3)
        # 输出 CPU times: user 318 ms, sys: 5.22 ms, total: 323 ms
        # Wall time: 337 ms
        # _theta: array([3.03182269, 3.93118623])

5.3 随机梯度下降法的封装和测试

    def fit_sgd(self, X_train, y_train, n_iters=5, t0=5, t1=50):
        """
        根据训练数据集X_train, y_train, 使用随机梯度下降法训练Linear Regression模型
        :param X_train:
        :param y_train:
        :param n_iters: 在随机梯度下降法中,n_iters代表所有的样本会被看几圈
        :param t0:
        :param t1:
        :return:
        """
        assert X_train.shape[0] == y_train.shape[0], \
            "the size of X_train must be equal to the size of y_train"
        assert n_iters >= 1

        def dJ_sgd(theta, X_b_i, y_i):
            """
            去X_b,y 中的随机一个元素进行导数公式的计算
            :param theta:
            :param X_b_i:
            :param y_i:
            :return:
            """
            return X_b_i * (X_b_i.dot(theta) - y_i) * 2.

        def sgd(X_b, y, initial_theta, n_iters, t0=5, t1=50):

            def learning_rate(t):
                """
                计算学习率,t1 为了减慢变化速度,t0为了增加随机性
                :param t: 第t次循环
                :return:
                """
                return t0 / (t + t1)

            theta = initial_theta
            m = len(X_b)

            for cur_iter in range(n_iters):
                # 对X_b进行一个乱序的排序
                indexes = np.random.permutation(m)
                X_b_new = X_b[indexes]
                y_new = y[indexes]

                # 对整个数据集看一遍
                for i in range(m):
                    gradient = dJ_sgd(theta, X_b_new[i], y_new[i])
                    theta = theta - learning_rate(cur_iter * m + i) * gradient

            return theta

        X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
        initial_theta = np.random.randn(X_b.shape[1])
        self._theta = sgd(X_b, y_train, initial_theta, n_iters, t0, t1)

        self.interception_ = self._theta[0]
        self.coef_ = self._theta[1:]

        return self

模拟数据进行测试

5.3-1

真实数据波士顿房价进行测试

5.3-2

5.4 使用Sklearn中的 随机梯度下降法

5.4-1

 

需要注意的是sklearn中的梯度下降法比我们自己的算法要复杂的多,性能和计算准确度上都比我们的要好,我们的算法只是用来演示过程,具体生产上的使用还是应该使用Sklearn提供的


6 梯度下降法 的调试

6.1 梯度下降法调试的原理

可能我们计算出梯度下降法的公式,并使用python编程实现,预测的过程中程序并没有报错,但是可能我们需要求的梯度的结果是错误的,这个时候需要怎么样去调试发现错误呢。

首先以二维坐标平面为例,一个点(O)的导数就是曲线在这个点的切线的斜率,在这个点两侧各取一个点(AB),那么AB两点对应的直线的斜率应该大体等于O的切线的斜率,并且这A和B的距离越近,那么两条直线的斜率就越接近
事实上,这也正是导数的定义,当函数y=f(x)的自变量x在一点x0上产生一个增量Δx时,函数输出值的增量Δy与自变量增量Δx的比值在Δx趋于0时的极限a如果存在,a即为在x0处的导数,记作f'(x0)或df(x0)/dx

6-1

 

扩展到多维维度则如下

6-2

梯度下降法调试的实现

np.random.seed(666)
X = np.random.normal(size=(1000,10))
X_b = np.hstack([np.ones((len(X),1)),X])

# 真实的θ值
true_theta = np.arange(1,12,dtype=float)
# np.random.normal(size=1000) 添加噪音
y = X_b.dot(true_theta) + np.random.normal(size=1000)
# 实现J(θ)函数
def J(theta,X_b,y):
    try:
        return np.sum((y-X_b.dot(theta))**2)/len(X_b)
    except:
        return float('inf')
    
# 实现数学推导出的dJ(θ)
def d_J_main(theta,X_b,y):
    return X_b.T.dot(X_b.dot(theta) - y) * 2. /len(X_b)

# 实现debug模式的dJ(θ)
def d_J_debug(theta,X_b,y,epsilon=0.01):
    res = np.empty(len(theta))
    for i in range(len(theta)):
        theta_1 = theta.copy()
        theta_1[i] += epsilon
        theta_2 = theta.copy()
        theta_2[i] -= epsilon
        res[i] = (J(theta_1,X_b,y)-J(theta_2,X_b,y))/(2*epsilon)
    return res
# 批量梯度下降法,d_J为求导函数,作为一个参数传入,用于切换求导策略
def gradient_descent(d_J,X_b, y, initial_theta, eta, n_iters=1e4, epsilon=1e-8):
            theta = initial_theta
            i_iters = 0

            while i_iters < n_iters:
                gradient = d_J(theta, X_b, y)
                last_theta = theta
                theta = theta - eta * gradient

                if abs(J(theta, X_b, y) - J(last_theta, X_b, y)) < epsilon:
                    break

                i_iters += 1

            return theta
X_b = np.hstack([np.ones((len(X), 1)), X])
initial_theta = np.zeros(X_b.shape[1])
eta = 0.01
# 使用d_J_debug调试模式求出theta
%time theta = gradient_descent(d_J_debug,X_b, y, initial_theta, eta)
print(theta)
# 使用数学解求出theta
%time theta = gradient_descent(d_J_main,X_b, y, initial_theta, eta)
print(theta)
# 输出结果
CPU times: user 531 ms, sys: 214 ms, total: 745 ms
Wall time: 613 ms
[ 0.94575233  1.98082712  3.06882065  3.94835863  4.97139932  5.9859077
  7.01077392  7.99250414  8.99151383  9.97525811 10.99758484]
CPU times: user 67 ms, sys: 27.6 ms, total: 94.6 ms
Wall time: 76.7 ms
[ 0.94575233  1.98082712  3.06882065  3.94835863  4.97139932  5.9859077
  7.01077392  7.99250414  8.99151383  9.97525811 10.99758484]

由此可以看出,我们的d_J_debug和d_J_main的结果是相近的,所以我们的d_J_main的数学推导是没问题的。
我们可以在真正的机器学习之前,先使用d_J_debug这种调试方式来验证一下我们的d_J_main的结果是否正确,然后再进行机器学习。

d_J_debug是通用的,可以放在任何求导的debug过程中,所以可以作为我们机器学习的工具箱来使用


7.梯度下降法的总结

7.1 小批量

  • 批量梯度下降法
  • 随机梯度下降法
    下面来看下二者的对比
维度 批量梯度下降法 随机梯度下降法
计算方式 每次对所有的样本看一遍才可以计算出梯度 每一次只需观察一个样本
速度
稳定性 高,一定可以先向损失函数下降的方式前进 低,每一次的方式不确定,甚至向反方向前进

综合二者的优缺点,有一种新的梯度下降法

7-1

 

小批量梯度下降法:即,我们每一次不看全部样本那么多,也不是只看一次样本那么少,每次只看k个样本
对于小批量梯度下降法,由多了一个超参数

def fit_lit_sgd(self, X_train, y_train, n_iters=5, t0=5, t1=50,k=10):
        """
        根据训练数据集X_train, y_train, 使用随机梯度下降法训练Linear Regression模型
        :param X_train:
        :param y_train:
        :param n_iters: 在随机梯度下降法中,n_iters代表所有的样本会被看几圈
        :param t0:
        :param t1:
        :param k: 小批量随机下降法的超参数k
        :return:
        """
        assert X_train.shape[0] == y_train.shape[0], \
            "the size of X_train must be equal to the size of y_train"
        assert n_iters >= 1

        def dJ_sgd(theta, X_b_k, y_k):
            """
            去X_b,y 中的随机选择k个元素进行导数公式的计算
            :param theta:
            :param X_b_i:
            :param y_i:
            :return:
            """
            return np.sum((X_b_k * (X_b_k.dot(theta) - y_k) ))* 2/len(X_b_k).

        def sgd(X_b, y, initial_theta, n_iters, t0=5, t1=50):

            def learning_rate(t):
                """
                计算学习率,t1 为了减慢变化速度,t0为了增加随机性
                :param t: 第t次循环
                :return:
                """
                return t0 / (t + t1)

            theta = initial_theta
            m = len(X_b)

            for cur_iter in range(n_iters):
                # 每次看k个元素
                i =0
                while i < m:
                    X_b_new = X_b[i:i+k]
                    y_new = y[i:i+k]
                    gradient = dJ_sgd(theta, X_b_new, y_new)
                    theta = theta - learning_rate(cur_iter * m + i+k) * gradient
                    i = i+k

            return theta

        X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
        initial_theta = np.random.randn(X_b.shape[1])
        self._theta = sgd(X_b, y_train, initial_theta, n_iters, t0, t1)

        self.interception_ = self._theta[0]
        self.coef_ = self._theta[1:]

        return self

7.2 随机

随机梯度下降法的优点

 

7.2-1

7.3 梯度上升法

7.3-1

7.3-2

7.3-3


 

版权声明:本文为博主原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。
本文链接:https://blog.csdn.net/qq_33182424/article/details/83108464

智能推荐

西门子HMI SMART 700 IE 设备概况以及WinCC flexible 2008常用配置小记-程序员宅基地

文章浏览阅读1.2w次,点赞3次,收藏32次。——参考自:SIMATIC HMI设备 Smart 700 IE、Smart 1000 IE 操作说明HMI 设备一旦探测到操作员控件被触摸就会立即返回一个反馈。该反馈是独立的,与 PLC 的通讯无关。 因此,其并不指示相关动作是否已真正执行。设计Smart Panel 700 IE是Smart Panel 700的升级版。连接HMI与PC连接组态PC与SMART PANEL有两种连接方式,一种是通过 RS485/422 接口连接,另一种是通过以太网接口连接,后者接线图如下图所示:HM_wincc flexible 2008

设备驱动模型(kobject、kset、ktype)_kobject kset ktype-程序员宅基地

文章浏览阅读570次。转自https://blog.csdn.net/guet_kite/article/details/78368928设备驱动模型概述Linux早期时候,一个驱动对应一个设备,也就对应一个硬件地址,那当有两个一样的设备的时候,就要写两个驱动,显然是不合理的。应该是从Linux2.5开始,就引入了device-bus-driver模型。其中设备驱动模型主要结构分为kset、kobject、kty..._kobject kset ktype

关于换行以及换行属性_nowarp也会换行-程序员宅基地

文章浏览阅读1.9k次。对于CSS的white-space属性,我想大部分人应该和我差不多,最常用的就是nowrap属性,最多用来做超长省略号显示的时候会用到【hiahiahia~】起因是这样的:产品doggie策划了一个元旦活动,活动主页最下边需要显示配置的活动规则,注意,是配置的活动规则,所以,免不了运营小妹要在后台配置一个活动规则,肯定不指望她们能配置html文本的呀,对吧,于是她们配置的是这样的: 1、封垫苏菲房间都是克拉夫; 2、对方萨芬的刷分放大; 3、粉打发打发打发这样色的,于..._nowarp也会换行

BRVAH万能适配器_andriod brvah适配器万能官网-程序员宅基地

文章浏览阅读353次。BRVAH一.简介二.使用一.简介BRVAH是一个强大的RecyclerAdapter框架(什么是RecyclerView?),它能节约开发者大量的开发时间,集成了大部分列表常用需求解决方案。二.使用在使用时,首先要项目的build.gradle导入allprojects { repositories { google() jcenter() ..._andriod brvah适配器万能官网

计算机文化基础作品ppt,计算机文化基础PPT课件-程序员宅基地

文章浏览阅读91次。计算机文化基础PPT课件2019-03-15计算机文化基础PPT课件 第1章http://wenku.baidu.com/view/e7ef8b6925c52cc58bd6be97.html计算机文化基础PPT课件 第2章 Windows 2000操作系统http://wenku.baidu.com/view/e85671f5f61fb7360b4c6594.html计算机文化基础PPT课件 ..._计算机文化基础ppt

获得迭代器最后一个元素_处理迭代器最后一个元素-程序员宅基地

文章浏览阅读3k次。来源python123获得迭代器最后一个元素问题尝试使用 * 迭代器展开运算,返回 range(0, 1000, 4) 的最后一个元素。print([x for x in range(0,1000,4)][-1])输出:996..._处理迭代器最后一个元素

随便推点

第七章 PX4-Pixhawk-Mavlink解析_px4 mavlink 波特率-程序员宅基地

文章浏览阅读5.2k次,点赞3次,收藏25次。第七章 PX4-Mavlink解析首先我们是还是来说一说mavlink吧。Mavlink协议是无人机的一种开源通信协议。可以理解就是按照一定的格式来发送数据。这一章节涉及到了消息的打包发送和接收解析。 首先我们还是找到入口函数然后回到脚本启动中找到mavlink的启动,这个找到应该不难吧,前面几章都有这个。这里有一个需要提一下,很多_px4 mavlink 波特率

Python 中RSA的用法 使用pyOpenssl 生成RSA密钥对, 使用rsa 加解密_import base64 import rsa from openssl.crypto impor-程序员宅基地

文章浏览阅读1.6k次,点赞4次,收藏8次。1. pyOpenSSL 生成RSA密钥对, 效率比较高2. rsa 加解密方法简单3. 代码import rsaimport base64from OpenSSL.crypto import PKeyfrom OpenSSL.crypto import TYPE_RSA, FILETYPE_PEM, FILETYPE_ASN1from OpenSSL.crypto import dum..._import base64 import rsa from openssl.crypto import pkey from openssl.crypto

Android 架构设计(四):组件化?_android 组件化 去除相关组件-程序员宅基地

文章浏览阅读1.8k次,点赞2次,收藏3次。同系列传送门Android 架构设计(一):设计模式分析_赵星海的博客-程序员宅基地Android 架构设计(二):分包和文件结构_赵星海的博客-程序员宅基地_android 分包结构Android 架构设计(三):技术选型_赵星海的博客-程序员宅基地关于组件化,我这边分三步与大家分享:1定义,2需求,3优劣,4改造步骤(含框架推荐);1、组件化的定义:各个业务模块可单独运行,模块相互联系只可以使用唯一的入口。如图:2、当前项目是否需要采用组件化?首先看项目大小,.._android 组件化 去除相关组件

通过设置偏移 添加RecyclerView分隔线_rv_list.additemdecoration(new recyclerviewdivider距-程序员宅基地

文章浏览阅读701次。添加RecyclerView分隔线_rv_list.additemdecoration(new recyclerviewdivider距离左侧

深入理解计算机系统--计算机系统漫游_深入理解计算机系统 jeancheng-程序员宅基地

文章浏览阅读282次。第一章 计算机系统漫游 计算机系统是由硬件和系统软件组成的。所有计算机系统都是由相似的硬件和软件组成,它们又执行着相似的功能。 以hello程序为例。 1.1信息就是位+上下文 hello程序的生命是从源程序(源文件)开始的。源程序是程序员编写的,hello.c。源程序是 0和1 的比特位,8个一组。ASCII标准来表示文本字符。 这样的文件称为文本文件,所有其他_深入理解计算机系统 jeancheng

python读取json字符串_json数据处理:读取文件中的json字符串,转为python字典-程序员宅基地

文章浏览阅读1.4k次。方法1:读取文件中的json字符串,再用json.loads转为python字典import jsonstr_file = ‘./960x540/config.json‘with open(str_file, ‘r‘) as f:print("Load str file from {}".format(str_file))str1 = f.read()r = json.loads(str1)pri..._python 提取json元素 获取两个字段的值组成字典