精华内容
下载资源
问答
  • 用C++实现了梯度下降求多元函数极值的算法,有可能会陷入局部最优解。
  • 梯度下降求极值

    千次阅读 2020-04-05 09:47:53
    梯度下降法的基本原理 梯度下降是迭代法的一种,可以用于求解最小二乘问题(线性和非线性都可以)。在求解机器学习算法的模型参数,即无约束优化问题时,梯度下降(Gradient Descent)是最常采用的方法之一,另一种常用...

    梯度下降法的基本原理

    梯度下降是迭代法的一种,可以用于求解最小二乘问题(线性和非线性都可以)。在求解机器学习算法的模型参数,即无约束优化问题时,梯度下降(Gradient Descent)是最常采用的方法之一,另一种常用的方法是最小二乘法。在求解损失函数的最小值时,可以通过梯度下降法来一步步的迭代求解,得到最小化的损失函数和模型参数值。反过来,如果我们需要求解损失函数的最大值,这时就需要用梯度上升法来迭代了。在机器学习中,基于基本的梯度下降法发展了两种梯度下降方法,分别为随机梯度下降法和批量梯度下降法。顾名思义,梯度下降法的计算过程就是沿梯度下降的方向求解极小值(也可以沿梯度上升方向求解极大值)。其迭代公式为,
    在这里插入图片描述
    其中代表梯负方向,表示梯度方向上的搜索步长。梯度方向我们可以通过对函数求导得到,步长的确定比较麻烦,太大了的话可能会发散,太小收敛速度又太慢。一般确定步长的方法是由线性搜索算法来确定,即把下一个点的坐标看做是ak+1的函数,然后求满足f(ak+1)的最小值的ak+1即可。因为一般情况下,梯度向量为0的话说明是到了一个极值点,此时梯度的幅值也为0.而采用梯度下降算法进行最优化求解时,算法迭代的终止条件是梯度向量的幅值接近0即可,可以设置个非常小的常数阈值。

    牛顿法

    牛顿迭代法是以微分为基础的,微分就是用直线来代替曲线,由于曲线不规则,那么我们来研究直线代替曲线后,剩下的差值是不是高阶无穷小,如果是高阶无穷小,那么这个差值就可以扔到不管了,只用直线就可以了,这就是微分的意义。
    牛顿法是牛顿在17世纪提出的一种求解方程f(x)=0.多数方程不存在求根公式,从而求精确根非常困难,甚至不可能,从而寻找方程的近似根就显得特别重要。
    牛顿迭代法是取x0之后,在这个基础上,找到比x0更接近的方程的跟,一步一步迭代,从而找到更接近方程根的近似跟。方法使用函数f(x)的泰勒级数的前面几项来寻找方程f(x) = 0的根。牛顿迭代法是求方程根的重要方法之一,其最大优点是在方程f(x) = 0的单根附近具有平方收敛,而且该法还可以用来求方程的重根、复根。另外该方法广泛用于计算机编程中。

    设r是f(x)=0的根,选取x0作为r初始近似值,过点(x0,f(x0))做曲线y=f(x)的切线L,L的方程为y=f(x0)+f’(x0)(x-x0),求出L与x轴交点的横坐标 x1=x0-f(x0)/f’(x0),称x1为r的一次近似值,过点(x1,f(x1))做曲线y=f(x)的切线,并求该切线与x轴的横坐标 x2=x1-f(x1)/f’(x1)称x2为r的二次近似值,重复以上过程,得r的近似值序列{Xn},其中Xn+1=Xn-f(Xn)/f’(Xn),称为r的n+1次近似值。上式称为牛顿迭代公式。

    梯度下降法求极值

    在这里插入图片描述

    在excel中求极值

    在这里插入图片描述
    在这里插入图片描述

    利用python代码求极值

    import math
    
    #使用梯度下降法求函数的最小值
    # f = x1^2+2x2^2-4x1-2x1x2  初始点为(1,1)
    
    
    
    #设计函数
    
    def function_one(x1,x2):   # 函数的输入 x,y
        f =  x1*x1+2*x2*x2-4*x1-2*x1*x2 # 算出f 的值
        dx = 2*x1-4-2*x2                          # 算出 一阶x导数的值
        dy = 4*x2-2*x1                      # 算出 一阶y导数的值
        
        return f, dx, dy                           # 返回三个值
    
    def main():
        x = 1 #初始点
        y = 1 #初始点
        
        a = 0.3  # 设置步长
        
        error = 0.000001  # 设置误差函数
        
        f_old, dx, dy = function_one(x,y) # 算出初始的值
        
        for i in range(100):    # 开始循环 梯度下降 设置循环的次数
            x -= a*dx
            y -= a*dy
            f_result, dx, dy = function_one(x,y)  # 获得第一次计算的初始值
            
            if abs(f_result - f_old) < error:       # 下降的结果绝对值与误差进行比较 如果再允许范围内则终止
                f_final = f_result
                print("极值点为: %.5f ,最小值为: %.5f, 循环次数: %d, 误差:%8f" % (x,f_final, i , abs(f_result - f_old)))
                break
            print("第 %d 次循环, 函数值为 %f" % (i, f_result))
            f_old = f_result
            
            
    if __name__ == '__main__':
        main()
    
    
    第 0 次循环, 函数值为 -5.400000
    第 1 次循环, 函数值为 -6.576000
    第 2 次循环, 函数值为 -7.193280
    第 3 次循环, 函数值为 -7.533504
    第 4 次循环, 函数值为 -7.727005
    第 5 次循环, 函数值为 -7.839158
    第 6 次循环, 函数值为 -7.904877
    第 7 次循环, 函数值为 -7.943626
    第 8 次循环, 函数值为 -7.966552
    第 9 次循环, 函数值为 -7.980142
    第 10 次循环, 函数值为 -7.988206
    第 11 次循环, 函数值为 -7.992994
    第 12 次循环, 函数值为 -7.995838
    第 13 次循环, 函数值为 -7.997527
    第 14 次循环, 函数值为 -7.998531
    第 15 次循环, 函数值为 -7.999127
    第 16 次循环, 函数值为 -7.999481
    第 17 次循环, 函数值为 -7.999692
    第 18 次循环, 函数值为 -7.999817
    第 19 次循环, 函数值为 -7.999891
    第 20 次循环, 函数值为 -7.999935
    第 21 次循环, 函数值为 -7.999962
    第 22 次循环, 函数值为 -7.999977
    第 23 次循环, 函数值为 -7.999986
    第 24 次循环, 函数值为 -7.999992
    第 25 次循环, 函数值为 -7.999995
    第 26 次循环, 函数值为 -7.999997
    第 27 次循环, 函数值为 -7.999998
    极值点为: 3.99862 ,最小值为: -8.00000, 循环次数: 28, 误差:0.000001
    

    由此可知,该函数的极值点为4.000,最小值为-8.000

    店铺多元回归

    先求方程和偏导

    #求第一组数据的
    from sympy import *
    a1,a2,b=symbols('a1 a2 b')
    S=(469-10*a1-80*a2-b)*(469-10*a1-80*a2-b)
    print(diff(S,a1))
    print(diff(S,a2))
    print(diff(S,b))
    
    200*a1 + 1600*a2 + 20*b - 9380
    1600*a1 + 12800*a2 + 160*b - 75040
    20*a1 + 160*a2 + 2*b - 938
    
    ##求第2组数据的
    from sympy import *
    a1,a2,b=symbols('a1 a2 b')
    S=(366-8*a1-0*a2-b)*(366-8*a1-0*a2-b)
    print(diff(S,a1))
    print(diff(S,a2))
    print(diff(S,b))
    
    128*a1 + 16*b - 5856
    0
    16*a1 + 2*b - 732
    
    ##求第3组数据的
    from sympy import *
    a1,a2,b=symbols('a1 a2 b')
    S=(371-8*a1-200*a2-b)*(371-8*a1-200*a2-b)
    print(diff(S,a1))
    print(diff(S,a2))
    print(diff(S,b))
    
    128*a1 + 3200*a2 + 16*b - 5936
    3200*a1 + 80000*a2 + 400*b - 148400
    16*a1 + 400*a2 + 2*b - 742
    
    ##求第4组数据的
    from sympy import *
    a1,a2,b=symbols('a1 a2 b')
    S=(208-5*a1-200*a2-b)*(208-5*a1-200*a2-b)
    print(diff(S,a1))
    print(diff(S,a2))
    print(diff(S,b))
    
    50*a1 + 2000*a2 + 10*b - 2080
    2000*a1 + 80000*a2 + 400*b - 83200
    10*a1 + 400*a2 + 2*b - 416
    
    ##求第5组数据的
    from sympy import *
    a1,a2,b=symbols('a1 a2 b')
    S=(246-7*a1-300*a2-b)*(246-7*a1-300*a2-b)
    print(diff(S,a1))
    print(diff(S,a2))
    print(diff(S,b))
    
    98*a1 + 4200*a2 + 14*b - 3444
    4200*a1 + 180000*a2 + 600*b - 147600
    14*a1 + 600*a2 + 2*b - 492
    
    ##求第6组数据的
    from sympy import *
    a1,a2,b=symbols('a1 a2 b')
    S=(297-8*a1-230*a2-b)*(297-8*a1-230*a2-b)
    print(diff(S,a1))
    print(diff(S,a2))
    print(diff(S,b))
    
    128*a1 + 3680*a2 + 16*b - 4752
    3680*a1 + 105800*a2 + 460*b - 136620
    16*a1 + 460*a2 + 2*b - 594
    
    ##求第7组数据的
    from sympy import *
    a1,a2,b=symbols('a1 a2 b')
    S=(363-7*a1-40*a2-b)*(363-7*a1-40*a2-b)
    print(diff(S,a1))
    print(diff(S,a2))
    print(diff(S,b))
    
    98*a1 + 560*a2 + 14*b - 5082
    560*a1 + 3200*a2 + 80*b - 29040
    14*a1 + 80*a2 + 2*b - 726
    
    ##求第8组数据的
    from sympy import *
    a1,a2,b=symbols('a1 a2 b')
    S=(436-9*a1-0*a2-b)*(436-9*a1-0*a2-b)
    print(diff(S,a1))
    print(diff(S,a2))
    print(diff(S,b))
    
    162*a1 + 18*b - 7848
    0
    18*a1 + 2*b - 872
    
    ##求第9组数据的
    from sympy import *
    a1,a2,b=symbols('a1 a2 b')
    S=(198-6*a1-330*a2-b)*(198-6*a1-330*a2-b)
    print(diff(S,a1))
    print(diff(S,a2))
    print(diff(S,b))
    
    72*a1 + 3960*a2 + 12*b - 2376
    3960*a1 + 217800*a2 + 660*b - 130680
    12*a1 + 660*a2 + 2*b - 396
    
    ##求第10组数据的
    from sympy import *
    a1,a2,b=symbols('a1 a2 b')
    S=(364-9*a1-180*a2-b)*(364-9*a1-180*a2-b)
    print(diff(S,a1))
    print(diff(S,a2))
    print(diff(S,b))
    
    162*a1 + 3240*a2 + 18*b - 6552
    3240*a1 + 64800*a2 + 360*b - 131040
    18*a1 + 360*a2 + 2*b - 728
    

    把结果加起来就是各自的偏导

    a1偏导=1226a1+22440a2+154b-53306
    a2偏导=22440a1+744400a2+3120b-881620
    b偏导=154a1+3120a2+20b-6636
    S=(469-10a1-80a2-b)(469-10a1-80a2-b)+(371-8a1-200a2-b)(371-8a1-200a2-b)+(371-8a1-200a2-b)(371-8a1-200a2-b)+(208-5a1-200a2-b)(208-5a1-200a2-b)+(246-7a1-300a2-b)(246-7a1-300a2-b)+(297-8a1-230a2-b)(297-8a1-230a2-b)+(363-7a1-40a2-b)(363-7a1-40a2-b)+(436-9a1-0a2-b)(436-9a1-0a2-b)+(198-6a1-330a2-b)(198-6a1-330a2-b)+(364-9a1-180a2-b)(364-9a1-180a2-b)

    在excel中求:

    在这里插入图片描述
    在这里插入图片描述

    由此可知:当a1=44.944,a2=-0.4475,b=69.9925时Se有极小值。
    所以y=44.944x1-0.4475x2+69.9925
    漫画书中所得结果y=41.5x1-0.3x2+65.3
    两者有一定的误差,不过用梯度下降法得到的结果没有用最小二乘法得到的结果准确。

    在Python中求

    import math
    
    #使用梯度下降法求函数的最小值
    
    #设计函数
    
    def function_one(x1,x2,b):   # 函数的输入 x,y
        f =  (469-10*x1-80*x2-b)*(469-10*x1-80*x2-b)+(366-8*x1-0*x2-b)*(366-8*x1-0*x2-b)+(371-8*x1-200*x2-b)*(371-8*x1-200*x2-b)+(208-5*x1-200*x2-b)*(208-5*x1-200*x2-b)+(246-7*x1-300*x2-b)*(246-7*x1-300*x2-b)+(297-8*x1-230*x2-b)*(297-8*x1-230*x2-b)+(363-7*x1-40*x2-b)*(363-7*x1-40*x2-b)+(436-9*x1-0*x2-b)*(436-9*x1-0*x2-b)+(198-6*x1-330*x2-b)*(198-6*x1-330*x2-b)+(364-9*x1-180*x2-b)*(364-9*x1-180*x2-b) # 算出f 的值
        dx1 = 1226*x1+22440*x2+154*b-53306                          # 算出 一阶x1导数的值
        dx2 = 22440*x1+744400*x2+3120*b-881620                     # 算出 一阶x2导数的值
        db = 154*x1+3120*x2+20*b-6636
        return f, dx1, dx2 ,db                          # 返回4个值
    
    def main():
        x1 = 45 #初始点
        x2 = 1 #初始点
        b = 70
        
        a = 0.0000001  # 设置步长
        
        error = 0.8  # 设置误差函数
        
        f_old, dx1, dx2,db = function_one(x1,x2,b) # 算出初始的值
        
        for i in range(100):    # 开始循环 梯度下降 设置循环的次数
            x1 -= a*dx1
            x2 -= a*dx2
            b -=b*db
            f_result, dx1, dx2,db = function_one(x1,x2,b)  # 获得第一次计算的初始值
            if abs(f_result-f_old)<error:
                print("a1=",x1)
                print("a2=",x2)
                print("b=",b)
            f_old=f_result
           
            
    if __name__ == '__main__':
        main()
    
    展开全文
  • 利用梯度下降函数极小值,或稍作修改用梯度上升法极大值,附带测试函数
  • 梯度与向量与梯度下降求极值

    千次阅读 2013-01-17 22:33:07
    梯度的模是在某个自变量变化为1时,应变量变化的量,如果是z=f(x,y),对x偏导,对y偏导,则梯度为:。例如的导数是2x,那么其意思是对应于应变量的变化情况,例如,对应x=0,对应的导数为0,联想图像,这点是...
           梯度是一个向量,一个方向,一维情况下是一维向量,二维情况下是二维向量。梯度的模是在某个自变量变化为1时,应变量变化的量,如果是求z=f(x,y),对x求偏导,对y求偏导,则梯度为:。例如的导数是2x,那么其意思是对应于应变量的变化情况,例如,对应x=0,对应的导数为0,联想图像,这点是没有增长的。如果x=1,那么导数对应2,即在x=1这点,当变化量为+1(相比这点增加1)时,应变量变化2。这里导数描述的是是变化率,导数为正,则变化1,对应增长,否则衰减。利用梯度下降法求函数的极小值是一种数值求解最优化问题的方法,对应利用梯度求最小值得问题,我理解的是假设原目标存存在最小的极值,梯度
    

    ,即对应某点的变化是方向中,指向变化最剧烈的方向,因为求极小值,则该方向对应导数为正,则反方向为极小值的方向,在penny梁斌的博客中,可以

    看到梯度下降法的具体介绍,还有BP算法解神经网络。

            当时看BP这篇的时候开始有点懵,理解上面的梯度就好理解了,在那个例子中,我们要求解的是各个神经元之间的参数w,在算导数是自然是对w求偏导。


    参考:BP算法浅谈 penny梁斌

    展开全文
  • 人工智能与机器学习——梯度下降函数极值一、原理介绍1. 梯度下降法的原理2. 梯度下降法求解过程3. 牛顿法原理4. 牛顿法方法说明二、用Excel完成函数极值的求解1. 求解函数题目2. 用excel计算① ∂z/∂x1 计算...

    一、原理介绍

    1. 梯度下降法的原理

    梯度下降是迭代法的一种,可以用于求解最小二乘问题(线性和非线性都可以)。在求解机器学习算法的模型参数,即无约束优化问题时,梯度下降(Gradient Descent)是最常采用的方法之一,另一种常用的方法是最小二乘法。在求解损失函数的最小值时,可以通过梯度下降法来一步步的迭代求解,得到最小化的损失函数和模型参数值。

    2. 梯度下降法求解过程

    顾名思义,梯度下降法的计算过程就是沿梯度下降的方向求解极小值(也可以沿梯度上升方向求解极大值)。

    其迭代公式为在这里插入图片描述, 其中在这里插入图片描述代表梯度负方向,
    在这里插入图片描述表示梯度方向上的搜索步长。梯度方向我们可以通过对函数求导得到,步长的确定比较麻烦,太大了的话可能会发散,太小收敛速度又太慢。一般确定步长的方法是由线性搜索算法来确定,即把下一个点的坐标看做是ak+1的函数,然后求满足f(ak+1)的最小值的ak+1即可。

    因为一般情况下,梯度向量为0的话说明是到了一个极值点,此时梯度的幅值也为0.而采用梯度下降算法进行最优化求解时,算法迭代的终止条件是梯度向量的幅值接近0即可,可以设置个非常小的常数阈值。

    3. 牛顿法原理

    牛顿法(英语:Newton’s method)又称为牛顿-拉弗森方法(英语:Newton-Raphson method),它是一种在实数域和复数域上近似求解方程的方法。方法使用函数在这里插入图片描述的泰勒级数的前面几项来寻找方程在这里插入图片描述的根。

    4. 牛顿法方法说明

    在这里插入图片描述

    二、用Excel完成函数极值的求解

    1. 求解函数题目

    用梯度法求二次函数的极小点及极小值,函数如下:二
    已知初始点在这里插入图片描述,迭代精度在这里插入图片描述

    2. 用excel计算

    在这里插入图片描述

    ① ∂z/∂x1 计算方法

    在这里插入图片描述

    ② ∂z/∂x2 计算方法

    在这里插入图片描述

    ③ ⊿x1 计算方法

    η*∂z/∂x1

    ④ ⊿x2 计算方法

    η*∂z/∂x2

    3. 计算结果

    在这里插入图片描述
    极小点(4,2) 极小值-8

    三、用 Python编程完成函数极值的求解

    1. 显示函数图像

    import numpy as np
    import matplotlib.pyplot as plt
    import matplotlib as mpl
    import math
    from mpl_toolkits.mplot3d import Axes3D
    import warnings
     
    # 解决中文显示问题
    mpl.rcParams['font.sans-serif'] = [u'SimHei']
    mpl.rcParams['axes.unicode_minus'] = False
    %matplotlib inline
     
    
    """
    二维原始图像
    1、构建一个函数为 y = x1^2 + 2 x2^2 - 4 x1 - 2 x1 x2 的图像。
    2、随机生成X1,X2点,根据X1,X2点生成Y点。
    3、画出图像。
    """
    def f2(x1,x2):
        return x1**2 + 2 * x2**2 - 4*x1 - 2 * x1 * x2 
     
    X1 = np.arange(-4,4,0.2)
    X2 = np.arange(-4,4,0.2)
    X1, X2 = np.meshgrid(X1, X2) # 生成xv、yv,将X1、X2变成n*m的矩阵,方便后面绘图
    Y = np.array(list(map(lambda t : f2(t[0],t[1]),zip(X1.flatten(),X2.flatten()))))
    Y.shape = X1.shape # 1600的Y图还原成原来的(40,40)
     
    %matplotlib inline
    #作图
    fig = plt.figure(facecolor='w')
    ax = Axes3D(fig)
    ax.plot_surface(X1,X2,Y,rstride=1,cstride=1,cmap=plt.cm.jet)
    ax.set_title(u'$ y = x1^2 + 2 x2^2 - 4 x1 - 2 x1 x2 $')
    plt.show()
    

    运行结果:
    在这里插入图片描述

    2. 求函数极值

    """
    对当前二维图像求最小点¶
    1、随机取一个点(x1,x2),设定α参数值  
    2、对这个点分别求关于x1、x2的偏导数,x1 =x1 - α*(dY/dx1),x2 =x2 - α*(dY/dx2)  
    3、重复第二补,设置 y的变化量 小于多少时 不再重复。
    """
    # 二维原始图像
    def f2(x, y):
        return x**2+2*(y)**2 - 4*(x)-2*(x)*(y)  
    ## 偏函数
    def hx1(x, y):
        return 2*x-4-2*y
    def hx2(x, y):
        return 4*y-2*x
    
    X1 = np.arange(-4,4,0.2)
    X2 = np.arange(-4,4,0.2)
    
    X1, X2 = np.meshgrid(X1, X2) # 生成xv、yv,将X1、X2变成n*m的矩阵,方便后面绘图
    Y = np.array(list(map(lambda t : f2(t[0],t[1]),zip(X1.flatten(),X2.flatten()))))
    Y.shape = X1.shape # 1600的Y图还原成原来的(40,40)
    
    
    x1 = 1
    x2 = 1
    alpha = 0.1
    #保存梯度下降经过的点
    GD_X1 = [x1]
    GD_X2 = [x2]
    GD_Y = [f2(x1,x2)]
    # 定义y的变化量和迭代次数
    y_change = f2(x1,x2)
    iter_num = 0
     
    while(y_change < 1e-10 and iter_num < 100) :
        tmp_x1 = x1 - alpha * hx1(x1,x2)
        tmp_x2 = x2 - alpha * hx2(x1,x2)
        tmp_y = f2(tmp_x1,tmp_x2)
        
        f_change = np.absolute(tmp_y - f2(x1,x2))
        x1 = tmp_x1
        x2 = tmp_x2
        GD_X1.append(x1)
        GD_X2.append(x2)
        GD_Y.append(tmp_y)
        iter_num += 1
    print(u"最终结果为:(%.5f, %.5f, %.5f)" % (x1, x2, f2(x1,x2)))
    print(u"迭代过程中X的取值,迭代次数:%d" % iter_num)
    print(GD_X1)
     
    # 作图
    fig = plt.figure(facecolor='w',figsize=(6,4))
    ax = Axes3D(fig)
    ax.plot_surface(X1,X2,Y,rstride=1,cstride=1,cmap=plt.cm.jet)
    ax.plot(GD_X1,GD_X2,GD_Y,'ko-')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.set_title(u'$ y = x1^2+2(x2)^2 - 4(x1)-2(x1) (x2) $')
     
    ax.set_title(u'函数;\n学习率:%.3f; 最终解:(%.3f, %.3f, %.3f);迭代次数:%d' % (alpha, x1, x2, f2(x1,x2), iter_num))
    plt.show()
    

    运行结果:
    在这里插入图片描述
    在这里插入图片描述

    四、梯度下降法求店铺面积与营业额问题

    1. excel表格

    在这里插入图片描述

    2. python代码

    import numpy as np
    import matplotlib.pyplot as plt
    import matplotlib as mpl
    from matplotlib import pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D
    
    # 解决中文显示问题
    mpl.rcParams['font.sans-serif'] = [u'SimHei']
    mpl.rcParams['axes.unicode_minus'] = False
    %matplotlib inline
    
    data=np.genfromtxt('E:\面积距离车站数据.csv',delimiter=',')
    x_data=data[:,:-1]
    y_data=data[:,2]
    #定义学习率、斜率、截据
    #设方程为y=theta1x1+theta2x2+theta0
    lr=0.00001
    theta0=0
    theta1=0
    theta2=0
    #定义最大迭代次数,因为梯度下降法是在不断迭代更新k与b
    epochs=10000
    #定义最小二乘法函数-损失函数(代价函数)
    def compute_error(theta0,theta1,theta2,x_data,y_data):
        totalerror=0
        for i in range(0,len(x_data)):#定义一共有多少样本点
            totalerror=totalerror+(y_data[i]-(theta1*x_data[i,0]+theta2*x_data[i,1]+theta0))**2
        return totalerror/float(len(x_data))/2
    #梯度下降算法求解参数
    def gradient_descent_runner(x_data,y_data,theta0,theta1,theta2,lr,epochs):
        m=len(x_data)
        for i in range(epochs):
            theta0_grad=0
            theta1_grad=0
            theta2_grad=0
            for j in range(0,m):
                theta0_grad-=(1/m)*(-(theta1*x_data[j,0]+theta2*x_data[j,1]+theta2)+y_data[j])
                theta1_grad-=(1/m)*x_data[j,0]*(-(theta1*x_data[j,0]+theta2*x_data[j,1]+theta0)+y_data[j])
                theta2_grad-=(1/m)*x_data[j,1]*(-(theta1*x_data[j,0]+theta2*x_data[j,1]+theta0)+y_data[j])
            theta0=theta0-lr*theta0_grad
            theta1=theta1-lr*theta1_grad
            theta2=theta2-lr*theta2_grad
        return theta0,theta1,theta2
    #进行迭代求解
    theta0,theta1,theta2=gradient_descent_runner(x_data,y_data,theta0,theta1,theta2,lr,epochs)
    print("多元线性回归方程为:y=",theta1,"X1+",theta2,"X2+",theta0)
    ax=plt.figure().add_subplot(111,projection='3d')
    ax.scatter(x_data[:,0],x_data[:,1],y_data,c='r',marker='o')
    x0=x_data[:,0]
    x1=x_data[:,1]
    #生成网格矩阵
    x0,x1=np.meshgrid(x0,x1)
    z=theta0+theta1*x0+theta2*x1
    #画3d图
    ax.plot_surface(x0,x1,z)
    ax.set_xlabel('店铺的面积',fontsize=20)
    ax.set_ylabel('距离最近的车站',fontsize=20)
    ax.set_zlabel("月营业额",fontsize=20)
    plt.show()
    

    3. 运行结果

    在这里插入图片描述

    五、最小二乘法求店铺面积与营业额问题

    1. excel表格

    在这里插入图片描述

    2. 求解二元线性回归方程

    from sklearn import linear_model        #表示,可以调用sklearn中的linear_model模块进行线性回归。
    import pandas as pd
    import seaborn as sns  
    import matplotlib.pyplot as plt
    %matplotlib inline
    
    data=pd.read_excel('E:\营业额数据.xlsx')
    
    X=data[['X1','X2']]   
    y=data['Y']    
    y=data.Y   
    
    model = linear_model.LinearRegression()
    model.fit(X,y)  
    

    在这里插入图片描述

    print("X1=",model.coef_[0])
    print("X2=",model.coef_[1])
    print("截距为=",model.intercept_)
    

    在这里插入图片描述

    print("多元线性回归方程为:y=",model.coef_[0],"X1",model.coef_[1],"X2+",model.intercept_)
    
    sns.pairplot(data, x_vars=['X1','X2'], y_vars='Y', height=3, aspect=0.8, kind='reg')  
    plt.show()  
    

    在这里插入图片描述

    六、比较最小二乘法和梯度下降法

    1. 梯度下降法

    在这里插入图片描述

    2. 最小二乘法

    在这里插入图片描述
    在这里插入图片描述

    七、参考文献

    《梯度下降》.天眼侍者.搜狗百科.2019-01-18
    链接: 梯度下降.

    《牛顿法》.≈风过无痕≈.搜狗百科.2018-10-25
    链接: 牛顿法.

    展开全文
  • MATLAB梯度下降多元函数的极值以及极值点; 程序+文档
  • 梯度下降梯度下降法的基本思想可以类比为一个下山的过程。 假设这样一个场景:一个人被困在山上,需要从山上下来(找到山的最低点,也就是山谷)。但此时山上的浓雾很大,导致可视度很低;因此,下山的路径就无法...
  • C++,二元函数,梯度下降法,求极值

    千次阅读 2018-05-14 16:17:24
    #include&lt;iostream&gt;#include &lt;iostream&gt;#include&lt;cmath&gt;#include &lt;random&gt; #include&lt;ctime&gt;using namespace std;... std::ra...
    #include<iostream>
    #include <iostream>
    #include<cmath>
    #include <random>
    #include<ctime>
    using namespace std;

    double normalf(double E, double sigema2)
    {
        double y;
        std::random_device rd{};
        std::mt19937 gen(rd());
        std::normal_distribution<> d{ E,sqrt(sigema2) }; //均值, 标准差
        //std::normal_distribution<> d{ E,sigema2};
        y = d(gen);
        return y;
    }

    //const double M_SQRT1_2 = sqrt(0.5);
    //double normalCFD(double value)  
    //{
    //    return 0.5 * erfc(-value * M_SQRT1_2);  // fai(x)
    //}

    const double M_SQRT1_2 = sqrt(0.5);
    double normpdpdf(double value)
    {
        return 0.5 * erfc(-value * M_SQRT1_2);
    }

    double qfunc(double value)
    {
        return 0.5 * erfc(value * M_SQRT1_2);  //  Q(x)     only  "-"
    }

    double f(double p, double q);
    double gp(double p, double q);
    double gq(double p, double q);
    double gd(double x0, double y0, double xx, double yx, double s);

    double gd(double x0, double y0, double xx, double yx, double s)
    {
        double x = x0;
        double y = y0; double Ipq;
        for (int i = 0; i != 40; ++i)
        {
            double grad1 = -1 * gp(x,y);
            double grad2 = -1 * gq(x, y);
            x += grad1*s; y += grad2*s;
            Ipq = f(x,y);//<< setprecision(15)
            cout << "EPOCH " << i << " grad1 = "  << grad1 << " x = " << x << "  y = " <<y <<"  Ipq = " <<Ipq << endl;
            if (abs(grad1)+ abs(grad2) < 1e-6)
                break;
        }
        xx = x;  yx = y;//x* ,y*
        return Ipq;
    }
    int main()
    {
        double xx=0; double yx=0;
        double x0 = -1,y0=4, ak = 0.1;
        gd(x0,y0,xx,yx, ak);
        for (;;);
    }

    double f(double p, double q )
    {    double IX;
        /*double u1, u2, sigema21, sigema22;
        double p1, p2, p3, p11, p21, p31, p1d, p2d, p3d, p11d, p21d, p31d;
        p1 = 1 - qfunc((p - u1) / sigema21);                       p1d = normpdpdf((p - u1) / sigema21);
        p2 = qfunc((p - u1) / sigema21) - qfunc((q - u1) / sigema21);  p2d = normpdpdf((q - u1) / sigema21) - normpdpdf((p - u1) / sigema21);
        p3 = qfunc((q - u1) / sigema21);                         p3d = -normpdpdf((q - u1) / sigema21);
        p11 = qfunc((q - u2) / sigema22);                        p11d = -normpdpdf((q - u2) / sigema22);
        p21 = qfunc((p - u2) / sigema22) - qfunc((q - u2) / sigema22); p21d = normpdpdf((q - u2) / sigema22) - normpdpdf((p - u2) / sigema22);
        p31 = 1 - qfunc((p - u2) / sigema22);                      p31d = normpdpdf((p - u2) / sigema22);    
        IX = 1 / 2 * (p1*log(2 * p1 / (p1 + p31)) + p2*log(2 * p2 / (p2 + p21)) + p3*log(2 * p3 / (p3 + p11)) + p11*log(2 * p11 / (p11 + p3)) + p21*log(2 * p21 / (p21 + p2)) + p31*log(2 * p31 / (p31 + p1)));
        */
        IX = p*p + q*q ;
        return IX;
    }
    double gp(double p, double q)
    {
        double DIXp;
        /*double u1, u2, sigema21, sigema22;
        double p1, p2, p3, p11, p21, p31, p1d, p2d, p3d, p11d, p21d, p31d;
        p1 = 1 - qfunc((p - u1) / sigema21);                       p1d = normpdpdf((p - u1) / sigema21);
        p2 = qfunc((p - u1) / sigema21) - qfunc((q - u1) / sigema21);  p2d = normpdpdf((q - u1) / sigema21) - normpdpdf((p - u1) / sigema21);
        p3 = qfunc((q - u1) / sigema21);                         p3d = -normpdpdf((q - u1) / sigema21);
        p11 = qfunc((q - u2) / sigema22);                        p11d = -normpdpdf((q - u2) / sigema22);
        p21 = qfunc((p - u2) / sigema22) - qfunc((q - u2) / sigema22); p21d = normpdpdf((q - u2) / sigema22) - normpdpdf((p - u2) / sigema22);
        p31 = 1 - qfunc((p - u2) / sigema22);                      p31d = normpdpdf((p - u2) / sigema22);

        DIXp = 1 / 2 * (p1d*log(2 * p1 / (p1 + p31)) + p2d*log(2 * p2 / (p2 + p21)) + p21d*log(2 * p21 / (p21 + p2)) + p31d*log(2 * p31 / (p31 + p1)));
        */
        DIXp = 2 * p;
        return DIXp;
    }
    double gq(double p, double q)
    {
        double DIXq;
        /*double u1, u2, sigema21, sigema22;
        double p1, p2, p3, p11, p21, p31, p1d, p2d, p3d, p11d, p21d, p31d;
        p1 = 1 - qfunc((p - u1) / sigema21);                       p1d = normpdpdf((p - u1) / sigema21);
        p2 = qfunc((p - u1) / sigema21) - qfunc((q - u1) / sigema21);  p2d = normpdpdf((q - u1) / sigema21) - normpdpdf((p - u1) / sigema21);
        p3 = qfunc((q - u1) / sigema21);                         p3d = -normpdpdf((q - u1) / sigema21);
        p11 = qfunc((q - u2) / sigema22);                        p11d = -normpdpdf((q - u2) / sigema22);
        p21 = qfunc((p - u2) / sigema22) - qfunc((q - u2) / sigema22); p21d = normpdpdf((q - u2) / sigema22) - normpdpdf((p - u2) / sigema22);
        p31 = 1 - qfunc((p - u2) / sigema22);                      p31d = normpdpdf((p - u2) / sigema22);
        DIXq = 1 / 2 * (p2d*log(2 * p2 / (p2 + p21)) + p3d*log(2 * p3 / (p3 + p11)) + p11d*log(2 * p11 / (p11 + p3)) + p21d*log(2 * p21 / (p21 + p2)));
        */
        DIXq = 2 * q;
        return DIXq;
    }

    展开全文
  • 一直疑惑一个问题:对于最小二乘法,为什么不直接求导让导数为“0”,直接求极值呢? 因为实际情况有些是不可行的,比如有时候求解这样的方程非常复杂。。 所以有了梯度下降:详细请查看:...
  • 泛函极值 变分法 梯度下降流1、积分泛函 欧拉-拉格朗日方程2、积分泛函 梯度下降流 在上一篇文章基于区域的主动轮廓模型-Chan Vese模型中,我们通过构造能量泛函,从而巧妙的把图像中的目标分割问题转化为求解能量...
  • 梯度下降初识-取凸函数极值 梯度下降梯度下降是迭代法的一种,通过选择一个初始点,然后计算该点的导数,再通过导数和步长推进到下一个点,直到两个点之间的差值很小为止。 # 梯度下降 == 导数值下降 import ...

空空如也

空空如也

1 2 3 4 5 ... 15
收藏数 292
精华内容 116
关键字:

梯度下降求极值