python实现共轭梯度法
作者:Tomator01 发布时间:2021-10-06 22:39:01
标签:python,共轭梯度法
共轭梯度法是介于最速下降法与牛顿法之间的一个方法,它仅需利用一阶导数信息,但克服了最速下降法收敛慢的缺点,又避免了牛顿法需要存储和计算Hesse矩阵并求逆的缺点,共轭梯度法不仅是解决大型线性方程组最有用的方法之一,也是解大型非线性最优化最有效的算法之一。 在各种优化算法中,共轭梯度法是非常重要的一种。其优点是所需存储量小,具有步收敛性,稳定性高,而且不需要任何外来参数。
算法步骤:
import random
import numpy as np
import matplotlib.pyplot as plt
def goldsteinsearch(f,df,d,x,alpham,rho,t):
'''
线性搜索子函数
数f,导数df,当前迭代点x和当前搜索方向d,t试探系数>1,
'''
flag = 0
a = 0
b = alpham
fk = f(x)
gk = df(x)
phi0 = fk
dphi0 = np.dot(gk, d)
alpha=b*random.uniform(0,1)
while(flag==0):
newfk = f(x + alpha * d)
phi = newfk
# print(phi,phi0,rho,alpha ,dphi0)
if (phi - phi0 )<= (rho * alpha * dphi0):
if (phi - phi0) >= ((1 - rho) * alpha * dphi0):
flag = 1
else:
a = alpha
b = b
if (b < alpham):
alpha = (a + b) / 2
else:
alpha = t * alpha
else:
a = a
b = alpha
alpha = (a + b) / 2
return alpha
def Wolfesearch(f,df,d,x,alpham,rho,t):
'''
线性搜索子函数
数f,导数df,当前迭代点x和当前搜索方向d
σ∈(ρ,1)=0.75
'''
sigma=0.75
flag = 0
a = 0
b = alpham
fk = f(x)
gk = df(x)
phi0 = fk
dphi0 = np.dot(gk, d)
alpha=b*random.uniform(0,1)
while(flag==0):
newfk = f(x + alpha * d)
phi = newfk
# print(phi,phi0,rho,alpha ,dphi0)
if (phi - phi0 )<= (rho * alpha * dphi0):
# if abs(np.dot(df(x + alpha * d),d))<=-sigma*dphi0:
if (phi - phi0) >= ((1 - rho) * alpha * dphi0):
flag = 1
else:
a = alpha
b = b
if (b < alpham):
alpha = (a + b) / 2
else:
alpha = t * alpha
else:
a = a
b = alpha
alpha = (a + b) / 2
return alpha
def frcg(fun,gfun,x0):
# x0是初始点,fun和gfun分别是目标函数和梯度
# x,val分别是近似最优点和最优值,k是迭代次数
# dk是搜索方向,gk是梯度方向
# epsilon是预设精度,np.linalg.norm(gk)求取向量的二范数
maxk = 5000
rho = 0.6
sigma = 0.4
k = 0
epsilon = 1e-5
n = np.shape(x0)[0]
itern = 0
W = np.zeros((2, 20000))
f = open("共轭.txt", 'w')
while k < maxk:
W[:, k] = x0
gk = gfun(x0)
itern += 1
itern %= n
if itern == 1:
dk = -gk
else:
beta = 1.0 * np.dot(gk, gk) / np.dot(g0, g0)
dk = -gk + beta * d0
gd = np.dot(gk, dk)
if gd >= 0.0:
dk = -gk
if np.linalg.norm(gk) < epsilon:
break
alpha=goldsteinsearch(fun,gfun,dk,x0,1,0.1,2)
# alpha=Wolfesearch(fun,gfun,dk,x0,1,0.1,2)
x0+=alpha*dk
f.write(str(k)+' '+str(np.linalg.norm(gk))+"\n")
print(k,alpha)
g0 = gk
d0 = dk
k += 1
W = W[:, 0:k+1] # 记录迭代点
return [x0, fun(x0), k,W]
def fun(x):
return 100 * (x[1] - x[0] ** 2) ** 2 + (1 - x[0]) ** 2
def gfun(x):
return np.array([-400 * x[0] * (x[1] - x[0] ** 2) - 2 * (1 - x[0]), 200 * (x[1] - x[0] ** 2)])
if __name__=="__main__":
X1 = np.arange(-1.5, 1.5 + 0.05, 0.05)
X2 = np.arange(-3.5, 4 + 0.05, 0.05)
[x1, x2] = np.meshgrid(X1, X2)
f = 100 * (x2 - x1 ** 2) ** 2 + (1 - x1) ** 2 # 给定的函数
plt.contour(x1, x2, f, 20) # 画出函数的20条轮廓线
x0 = np.array([-1.2, 1])
x=frcg(fun,gfun,x0)
print(x[0],x[2])
# [1.00318532 1.00639618]
W=x[3]
# print(W[:, :])
plt.plot(W[0, :], W[1, :], 'g*-') # 画出迭代点收敛的轨迹
plt.show()
代码中求最优步长用得是goldsteinsearch方法,另外的Wolfesearch是试验的部分,在本段程序中不起作用。
迭代轨迹:
三种最优化方法的迭代次数对比:
最优化方法 | 最速下降法 | 共轭梯度法 | 牛顿法 |
迭代次数 | 1702 | 240 | 5 |
来源:https://blog.csdn.net/Big_Pai/article/details/88540147
0
投稿
猜你喜欢
- 一:数据源介绍本篇文章介绍的是使用python实现对葵花8号卫星数据进行自动下载。葵花8号卫星是日本的一颗静止轨道气象卫星,覆盖范围为60S
- 数据库是数据的集合,与数学的集合论有密不可分的关系。为提高查询速度,我们可以:对数据表添加索引,以加快搜索速度;通过编程技巧最大限度地利用索
- 连接数据库:mysql -u用户名 -p密码导入数据 source d:\create.sql用下面的语句就可以导出mysql中的数据了:m
- 从大的发展来看,网站就是一块试验田,一块在错误中成长、在错误中变强变大的试验田。这决定了互联网产品的成长路线,一定是一个反复修正和迭代的曲线
- 下拉菜单平常见到的都是用js来实现的,本文介绍的方法是使用纯CSS实现导航下拉菜单功能,代码符合标准,兼容性好且环保,制作下拉菜单的不错选择
- 步骤一:下载对应的CURL压缩包并在windows上配置好环境变量进入CURL官网下载对应的windows压缩包。地址:点击打开链接把下载好
- 如下所示:"""提取文档数超过10000的数据按照某个字段的值具有唯一性进行升序,按照@timestamp进行
- Telnet 连接方式#!/usr/bin/env python# coding=utf-8import timeimport telnet
- 随着信息技术的飞速发展,数据处理不仅在数量上要求越来越大,而且在质量上也要求越来越高。操作系统的稳定对数据库来说是十分紧要的,在数据库可操作
- 1、实现效果2、实现步骤模块导入import os,sys,timefrom PyQt5 import QtCore,QtWidgets,Q
- insert into(列名) select 列名 from 表名 where 条件 --不创建表,只复制表数据 select 列名 int
- <% If Err.Number <> 0 Th
- 摘要:本文介绍了字符与编码的发展过程,相关概念的正确理解。举例说明了一些实际应用中,编码的实现方法。然后,本文讲述了通常对字符与编码的几种误
- ObjectUtil组件其实就是单例模式的最好范例,声明调用各个组件的时候,用ObjectUtil调用,可以有效的防止调用过多的类而导致错误
- 一、安装cuda1、在英伟达官网下载最新版的cuda驱动https://developer.nvidia.com/zh-cn/cuda-do
- IE8主页http://www.microsoft.com/windows/products/winfamily/ie/ie8/defaul
- urllib的基本用法urllib库的基本组成利用最简单的urlopen方法爬取网页html利用Request方法构建headers模拟浏览
- 一、简单使用def TestLogBasic(): import logging  
- 本文实例讲述了Python使用迭代器捕获Generator返回值的方法。分享给大家供大家参考,具体如下:用for循环调用generator时
- 鼠标双击滚动屏幕,单击停止滚动,很多小说新闻网站都有这个很人性化的功能,阅读起小说、新闻来很方便,不用手动拉滚动条。js代码如下:<h