IT数码 购物 网址 头条 软件 日历 阅读 图书馆
TxT小说阅读器
↓语音阅读,小说下载,古典文学↓
图片批量下载器
↓批量下载图片,美女图库↓
图片自动播放器
↓图片自动播放器↓
一键清除垃圾
↓轻轻一点,清除系统垃圾↓
开发: C++知识库 Java知识库 JavaScript Python PHP知识库 人工智能 区块链 大数据 移动开发 嵌入式 开发工具 数据结构与算法 开发测试 游戏开发 网络协议 系统运维
教程: HTML教程 CSS教程 JavaScript教程 Go语言教程 JQuery教程 VUE教程 VUE3教程 Bootstrap教程 SQL数据库教程 C语言教程 C++教程 Java教程 Python教程 Python3教程 C#教程
数码: 电脑 笔记本 显卡 显示器 固态硬盘 硬盘 耳机 手机 iphone vivo oppo 小米 华为 单反 装机 图拉丁
 
   -> 人工智能 -> python实现共轭梯度算法 -> 正文阅读

[人工智能]python实现共轭梯度算法

一、共轭梯度算法简介

共轭梯度(Conjugate Gradient)方法是一种迭代算法,可用于求解无约束的优化问题,例如能量最小化。常见的优化算法还有梯度下降法,相比于梯度下降法,共轭梯度法具有收敛更快,以及算法稳定性更好的优点。
梯度下降法和共轭梯度法
从上图可以看出来,梯度下降法优化过程中函数是沿着梯度的反方向逐步优化,后一步优化的结果会对前一步的优化结果造成影响,收敛较慢。而共轭梯度方法每一步的优化方向与前一步的优化方向是共轭的,因此并不会对前一步的优化结果造成影响,同时优化过程中保证在每一个方向上函数优化到最小值,从而保证沿着这些共轭方向优化完成后,函数能够达到全局最小值点。具体解释可以参考这篇文章:https://zhuanlan.zhihu.com/p/338838078

二、实现共轭梯度方法的两块重要积木

1.共轭方向的确定

共轭梯度方法中的新的共轭方向仅由上一步的梯度新的梯度上一步的优化方向决定,即:
d ? k + 1 = ? g ? k + 1 + β k d ? k \vec{d}_{k+1}=-\vec{g}_{k+1}+\beta_k\vec{d}_k d k+1?=?g ?k+1?+βk?d k?
其中 β k \beta_k βk?的定义方式有很多种,常用的有FRPRP两种:
β k F R = g k + 1 T g k + 1 g k T g k = ∣ g ? k + 1 ∣ 2 ∣ g ? k ∣ 2 \beta_k^{FR}=\frac{g_{k+1}^Tg_{k+1}}{g_k^Tg_k}=\frac{|\vec{g}_{k+1}|^2}{|\vec{g}_k|^2} βkFR?=gkT?gk?gk+1T?gk+1??=g ?k?2g ?k+1?2?
β k P R P = g k + 1 T ( g k + 1 ? g k ) g k T g k = g ? k + 1 ? ( g ? k + 1 ? g ? k ) ∣ g ? k ∣ 2 \beta_k^{PRP}=\frac{g_{k+1}^T(g_{k+1}-g_k)}{g_k^Tg_k}=\frac{\vec{g}_{k+1}\cdot(\vec{g}_{k+1}-\vec{g}_k)}{|\vec{g}_k|^2} βkPRP?=gkT?gk?gk+1T?(gk+1??gk?)?=g ?k?2g ?k+1??(g ?k+1??g ?k?)?

2.方向优化步长的确定

当优化方向确定之后,需要利用线性搜索技术确定优化的步长 α m i n \alpha_{min} αmin?,即是寻找α>0使得
f ( x ? + α m i n d ? ) = m i n f ( x ? + α d ? ) f(\vec{x}+\alpha_{min}\vec{d})=minf(\vec{x}+\alpha\vec{d}) f(x +αmin?d )=minf(x +αd )
线搜索技术包含两大类,即精确线搜索技术和非精确线搜索技术,详见这篇文章:https://www.longzf.com/optimization/2/line_search/
精确线搜索技术包括牛顿法和二分法。由于之后的代码使用的线搜索技术是牛顿法,所以下面简单介绍一下牛顿法。
牛顿法被用于求解函数的极小极大值问题,函数在极值点处应有 f ′ ( α ) = 0 f'(\alpha)=0 f(α)=0,将函数展开到二阶泰勒展开之后,可以得到α的迭代公式:
α k + 1 = α k ? f ′ ( α k ) f ′ ′ ( α k ) \alpha_{k+1} = {\alpha}_k - \frac{f'(\alpha_k)}{f''(\alpha_k)} αk+1?=αk??f(αk?)f(αk?)?

note

由于数值计算过程中要用到差分的方法,所以这里简单列出以下使用一阶微分和二阶微分所使用的差分表达式:
f ′ ( x ) ≈ f ( x + δ ) ? f ( x ? δ ) 2 δ f'(x)\approx\frac{f(x+\delta )-f(x-\delta)}{2\delta} f(x)2δf(x+δ)?f(x?δ)?
f ′ ′ ( x ) ≈ f ( x + δ ) + f ( x ? δ ) ? 2 f ( x ) δ 2 f''(x)\approx\frac{f(x+\delta)+f(x-\delta)-2f(x)}{\delta^2} f(x)δ2f(x+δ)+f(x?δ)?2f(x)?
其中 δ \delta δ是小量。

三、共轭梯度算法优化过程

  1. 计算初始初始梯度值 g ? 0 \vec{g}_0 g ?0?和优化方向 d ? 0 \vec{d}_0 d 0?
    g ? 0 = ? ? f ( x ? 0 ) \vec{g}_0=\vec{\nabla}f(\vec{x}_0) g ?0?=? f(x 0?)
    d ? 0 = ? g ? 0 \vec{d}_0 = -\vec{g}_0 d 0?=?g ?0?

  2. 如果 g ? k < ? ? \vec{g}_k<\vec{\epsilon} g ?k?<? ,退出迭代过程,否则执行以下步骤

  3. 用线性搜索算法(牛顿法)求出使得 m i n f ( x ? k + α d ? k ) minf(\vec{x}_k+\alpha\vec{d}_k) minf(x k?+αd k?)的步长α,并更新 x ? k + 1 \vec{x}_{k+1} x k+1?:
    x ? k + 1 = x ? k + α d ? k \vec{x}_{k+1}=\vec{x}_k+\alpha\vec{d}_k x k+1?=x k?+αd k?

  4. 计算新的梯度
    g ? k + 1 = ? ? f ( x ? k + 1 ) \vec{g}_{k+1}=\vec{\nabla}f(\vec{x}_{k+1}) g ?k+1?=? f(x k+1?)

  5. 根据前面提到的FR公式或PRP公式计算新的组合系数

  6. 计算新的共轭方向:
    d ? k + 1 = ? g ? k + 1 + β k d ? k \vec{d}_{k+1}=-\vec{g}_{k+1}+\beta_k\vec{d}_k d k+1?=?g ?k+1?+βk?d k?

  7. 重复执行第2步

四、python实现共轭梯度算法

这里使用的测试函数形式是:
f ( x , y ) = ( 3 x ? 2 y ) 2 + ( x ? 1 ) 4 f(x, y)=(3x-2y)^2+(x-1)^4 f(x,y)=(3x?2y)2+(x?1)4
可以看到函数存在最小值点(1, 1.5),以下是实现的python代码。

1.FR-CG

############共轭梯度算法#############
import numpy as np
import matplotlib.pyplot as plt
def testFun(x, y):
    t = 3.0*x - 2.0*y
    t1 = x - 1.0
    z = np.power(t, 2) + np.power(t1, 4)
    return z

def gradTestFun(x, y):
    '''
    求函数的梯度
    '''
    delta_x = 1e-6      #x方向差分小量
    delta_y = 1e-6      #y方向差分小量
    grad_x = (testFun(x+delta_x, y)-testFun(x-delta_x, y))/(2.0*delta_x)
    grad_y = (testFun(x, y+delta_y)-testFun(x, y-delta_y))/(2.0*delta_y)
    grad_xy = np.array([grad_x, grad_y])
    return grad_xy

def getStepLengthByNewton(array_xy, array_d):
    '''
    采用牛顿法,精确线性搜索确定移动步长
    '''
    a0 = 1.0           #初始猜测值
    e0 = 1e-6          #退出搜索循环的条件
    delta_a = 1e-6     #对a作差分的小量
    while(1):
        new_a = array_xy + a0*array_d
        new_a_l = array_xy + (a0-delta_a)*array_d
        new_a_h = array_xy + (a0+delta_a)*array_d
        diff_a0 = (testFun(new_a_h[0], new_a_h[1]) - testFun(new_a_l[0], new_a_l[1]))/(2.0*delta_a)
        if np.abs(diff_a0) < e0:
            break
        ddiff_a0 = (testFun(new_a_h[0], new_a_h[1]) + testFun(new_a_l[0], new_a_l[1]) - 2.0*testFun(new_a[0], new_a[1]))/(delta_a*delta_a)
        a0 = a0 - diff_a0/ddiff_a0
    return a0

def plotResult(array_xy_history):
    x = np.linspace(-1.0, 4.0, 100)
    y = np.linspace(-4.0, 8.0, 100)
    X, Y = np.meshgrid(x, y)
    Z = testFun(X, Y)
    plt.figure(dpi=300)
    plt.xlim(-1.0, 4.0)
    plt.ylim(-4.0, 8.0)
    plt.xlabel("x")
    plt.ylabel("y")
    plt.contour(X, Y, Z, 40)
    plt.plot(array_xy_history[:,0], array_xy_history[:,1], marker='.', ms=10)
    xy_count = array_xy_history.shape[0]
    for i in range(xy_count):
        if i == xy_count-1:
            break
        dx = (array_xy_history[i+1][0] - array_xy_history[i][0])*0.6
        dy = (array_xy_history[i+1][1] - array_xy_history[i][1])*0.6
        plt.arrow(array_xy_history[i][0], array_xy_history[i][1], dx, dy, width=0.1)

def mainFRCG():
    '''
    使用CG算法优化,用FR公式计算组合系数
    '''
    ls_xy_history = []                           #存储初始坐标的迭代结果
    xy0 = np.array([4.0, -2.0])                   #初始点
    grad_xy = gradTestFun(xy0[0], xy0[1])
    d = -1.0*grad_xy                             #初始搜索方向
    e0 = 1e-6                                    #迭代退出条件
    xy = xy0
    while(1):
        ls_xy_history.append(xy)
        grad_xy = gradTestFun(xy[0], xy[1])
        tag_reach = np.abs(grad_xy) < e0
        if tag_reach.all():
            break
        step_length = getStepLengthByNewton(xy, d)
        xy_new = xy + step_length*d
        grad_xy_new = gradTestFun(xy_new[0], xy_new[1])
        b = np.dot(grad_xy_new, grad_xy_new)/np.dot(grad_xy, grad_xy)       #根据FR公式计算组合系数
        d = b*d - grad_xy_new
        xy = xy_new
    array_xy_history = np.array(ls_xy_history)
    plotResult(array_xy_history)
    return array_xy_history

以下是运行结果
在这里插入图片描述

最终得到的最小值的坐标为[ 1.00113664, 1.50170496]
运行时间:
在这里插入图片描述

2.PRP-CG

换用PRP公式计算只需要将上面求解组合系数部分的代码作少量修改即可,代码为:

def mainPRPCG():
    '''
    使用CG算法优化,用PRP公式计算组合系数
    '''
    ls_xy_history = []
    xy0 = np.array([4.0, -2.0])
    grad_xy = gradTestFun(xy0[0], xy0[1])
    d = -1.0*grad_xy
    e0 = 1e-6
    xy = xy0
    while(1):
        ls_xy_history.append(xy)
        grad_xy = gradTestFun(xy[0], xy[1])
        tag_reach = np.abs(grad_xy) < e0
        if tag_reach.all():
            break
        step_length = getStepLengthByNewton(xy, d)
        xy_new = xy + step_length*d
        grad_xy_new = gradTestFun(xy_new[0], xy_new[1])
        b = np.dot(grad_xy_new, (grad_xy_new - grad_xy))/np.dot(grad_xy, grad_xy)    #根据PRP公式计算组合系数
        d = b*d - grad_xy_new
        xy = xy_new
    array_xy_history = np.array(ls_xy_history)
    plotResult(array_xy_history)
    return array_xy_history

运行结果为:
在这里插入图片描述
得到的最小值的坐标是[ 1.00318321, 1.50477481]
运行时间:
在这里插入图片描述

3.GD

最后还一起写了利用梯度下降法优化的结果以作对比,以下是实现代码:

def mainGD():
    '''
    使用梯度下降法计优化函数
    '''
    ls_xy_history = []
    xy0 = np.array([4.0, -2.0])
    ls_xy_history.append(xy0)
    grad_xy = gradTestFun(xy0[0], xy0[1])
    alpha = 1e-3      
    e0 = 1e-3
    xy = xy0
    while(1):
        tag_reach = np.abs(grad_xy)<e0
        if tag_reach.all():
            break
        xy = xy - alpha*grad_xy
        grad_xy = gradTestFun(xy[0], xy[1])
        ls_xy_history.append(xy)
    array_xy_history = np.array(ls_xy_history)
    plotResult(array_xy_history)
    return array_xy_history

运行结果为:
在这里插入图片描述
得到的最小值坐标是[ 0.9185095 , 1.37763925]
运行时间是:
在这里插入图片描述

五、总结

对比可以发现,CG算法在计算的速度和准确度上都相较于GD算法有一定的优势。

  人工智能 最新文章
2022吴恩达机器学习课程——第二课(神经网
第十五章 规则学习
FixMatch: Simplifying Semi-Supervised Le
数据挖掘Java——Kmeans算法的实现
大脑皮层的分割方法
【翻译】GPT-3是如何工作的
论文笔记:TEACHTEXT: CrossModal Generaliz
python从零学(六)
详解Python 3.x 导入(import)
【答读者问27】backtrader不支持最新版本的
上一篇文章      下一篇文章      查看所有文章
加:2022-03-24 00:32:32  更:2022-03-24 00:36:13 
 
开发: C++知识库 Java知识库 JavaScript Python PHP知识库 人工智能 区块链 大数据 移动开发 嵌入式 开发工具 数据结构与算法 开发测试 游戏开发 网络协议 系统运维
教程: HTML教程 CSS教程 JavaScript教程 Go语言教程 JQuery教程 VUE教程 VUE3教程 Bootstrap教程 SQL数据库教程 C语言教程 C++教程 Java教程 Python教程 Python3教程 C#教程
数码: 电脑 笔记本 显卡 显示器 固态硬盘 硬盘 耳机 手机 iphone vivo oppo 小米 华为 单反 装机 图拉丁

360图书馆 购物 三丰科技 阅读网 日历 万年历 2025年1日历 -2025/1/9 1:22:05-

图片自动播放器
↓图片自动播放器↓
TxT小说阅读器
↓语音阅读,小说下载,古典文学↓
一键清除垃圾
↓轻轻一点,清除系统垃圾↓
图片批量下载器
↓批量下载图片,美女图库↓
  网站联系: qq:121756557 email:121756557@qq.com  IT数码