• 主要介绍了python实现的共轭梯度法，文中通过示例代码介绍的非常详细，对大家的学习或者工作具有一定的参考学习价值，需要的朋友们下面随着小编来一起学习学习吧
• CG共轭梯度算法python实现 # 775.5289919376373 500维度 def CG2(A, b, x, imax=500, epsilon=0.0000001): steps = np.asarray(x) i = 0 # 计算残差，也是负梯度方向 r = b - np.dot(A ,x)
• 资源包括两个算法的python实现，分别是共轭梯度法和BFGS法，利用Numpy与Sympy，实现了只需输入函数与初始点等基本条件，即可求解并输出迭代过程种各参数变化。
• from numpy import * import numpy import numpy as np def f(A,x,b): return 0.5*np.dot(np.dot(x,A),x)+np.dot(b,x) def g(A,x,b): return np.dot(A,x)+b # def minimize_cg(A,b, x0, jac=None,gtol=1e...

from numpy import *
import numpy
import numpy as np

def f(A,x,b):
return 0.5*np.dot(np.dot(x,A),x)+np.dot(b,x)
def  g(A,x,b):
return np.dot(A,x)+b
#
def minimize_cg(A,b, x0, jac=None,gtol=1e-5,maxiter=None,disp=False):
"""
Minimization of scalar function of one or more variables using the

Options
-------
disp : bool
Set to True to print convergence messages.
maxiter : int
Maximum number of iterations to perform.
gtol : float
Gradient norm must be less than gtol before successful
termination.
norm : float
Order of norm (Inf is max, -Inf is min).

"""
if maxiter is None:
maxiter = len(x0) * 200
gfk = np.dot(A, x0) + b
k = 0
xk = x0
warnflag = 0
pk = -gfk

gnorm = numpy.amax(numpy.abs(gfk))
while (gnorm > gtol) and (k < maxiter):
deltak = numpy.dot(gfk, gfk)
alpha_k = -np.dot(gfk, pk) / (np.dot(np.dot(pk,A.T ),pk))
xk = xk + alpha_k * pk
gfkp1=np.dot(A, xk) + b
beta_k = max(0, numpy.dot(gfkp1, gfkp1) / deltak)
pk = -gfkp1 + beta_k * pk
gfk = gfkp1
gnorm=numpy.amax(numpy.abs(gfk))//最大值作为范数
k += 1
if warnflag == 0:
# msg = _status_message['success']
if disp:
print('success')
print("         Current function value: " , xk )
print("         Iterations: %d" % k)
print("         Function evaluations: " , f(A,xk,b))
print("         actual value x*: ",-np.dot(np.linalg.inv(A),b))
print("         actual function value f(x*): ",f(A,-np.dot(np.linalg.inv(A),b),b))

if __name__ == '__main__':
x0 = [2, 1]
A=np.array([[4,1],[1,3]])
b=[-1,-2]
g0=g(A,x0,b)
minimize_cg(A,b,x0,g0,disp=True,maxiter=100)


展开全文
• 主要介绍了基于Python共轭梯度法与最速下降法之间的对比，具有很好的参考价值，希望对大家有所帮助。一起跟随小编过来看看吧
• matlab共轭梯度法求目标函数的最小极值-共轭梯度-王.rar 我是地球物理专业的一名学生，把自己实习的作业发上来大家分享下
• Python实现最速下降法、共轭梯度法和信赖域狗腿法源代码。可以直接运行，同时将迭代分析绘图。配有详细注释
• 共轭梯度法实现求解线性方程组。 可以用来求解一般的线性方程组方程，程序清晰易懂。
• 第十二课 共轭梯度法解方程 运行上一篇的程序得到的结果表明，最陡下降法与迄今为止所其它迭代方法比较，在迭代次数上没有竞争力。然而，如果从根本上改进与[A]相互“共轭”的下降向量。于是，引入了满足以下关系的...
第十二课 共轭梯度法解方程
运行上一篇的程序得到的结果表明，最陡下降法与迄今为止所其它迭代方法比较，在迭代次数上没有竞争力。然而，如果从根本上改进与[A]相互“共轭”的下降向量。于是，引入了满足以下关系的“下降向量” 按照我的个人理解，共轭梯度法不仅在试解上做修正，还在误差值上做改进，这样做的目的是为了提高迭代效率  最陡下降法将被改进为，
算例依然采用高斯赛德尔的例子 程序代码如下：分别为一个主程序和一个检查收敛的子程序checkit 主程序
#线性联立方程的共轭梯度法
import numpy as np
import math
import B
n=3
converged=np.array([False])
p=np.zeros((n,1))
r=np.zeros((n,1))
u=np.zeros((n,1))
xnew=np.zeros((n,1))
a=np.array([[16,4,8],[4,5,-4],[8,-4,22]],dtype=np.float)
b=np.array([[4],[2],[5]],dtype=np.float)
x=np.array([[1],[1],[1]],dtype=np.float)
tol=1.0e-5
limit=100
print('系数矩阵')
print(a[:])
print('右手边向量',b[:,0])
print('初始猜测值',x[:,0])
r[:]=b[:]-np.dot(a,x)
p[:]=r[:]
print('前几次迭代值')
iters=0
while(True):
iters=iters+1
u[:]=np.dot(a,p)
up=np.dot(np.transpose(r),r)
alpha=up/np.dot(np.transpose(p),u)
xnew[:]=x[:]+p[:]*alpha
r[:]=r[:]-u[:]*alpha
beta=np.dot(np.transpose(r),r)/up
p[:]=r[:]+p[:]*beta
if iters<5:
print(x[:,0])
B.checkit(xnew,x,tol,converged)
if converged==True or iters==limit:
break
x[:,0]=xnew[:,0]
print('到收敛需要迭代次数',iters)
print('解向量',x[:,0])

checkit

def checkit(loads,oldlds,tol,converged):
#检查前后两个量的收敛性
big=0.0
converged[:]=True
for i in range(1,neq+1):
for i in range(1,neq+1):
converged[:]=False

终端输出结果：  程序结果与计算结果一致，而且只需要4次迭代，可见共轭梯度法的高效性。
展开全文
• import numpy as np def CG(A,b): n = b.shape[0] xs = [] rs = [] ps = [] alphas = [] x0 = np.array([2,1]) # x0 = np.random.rand(n) xs.append(x0) r0 = b - np.dot(A,x0...
import numpy as np
def CG(A,b):
n = b.shape[0]
xs = []
rs = []
ps = []
alphas = []
x0 = np.array([2,1])
# x0 = np.random.rand(n)
xs.append(x0)

r0 = b - np.dot(A,x0)
rs.append(r0)

p0 = r0
ps.append(p0)

alpha0 = p0.dot(p0)/p0.dot(A).dot(p0)
alphas.append(alpha0)
print(rs)
print(alphas)
print(ps)
for i in range(n):
r = rs[i] - alphas[i] * A.dot(ps[i])
rs.append(r)
beta = np.dot(r,r)/(rs[i].dot(rs[i]))

alpha = ps[i].dot(rs[i])/(ps[i]).dot(A).dot(ps[i])
alphas.append(alpha)

x = xs[i] + alpha * ps[i]
xs.append(x)

p = r + beta * ps[i]
ps.append(p)

return xs

A = np.array([[4,1],[1,3]])
b = np.array([1,2])
c1 = np.linalg.inv(A).dot(b)
print("the math sol is",c1)
c2 = CG(A,b)
print("the numerical sol is",c2)


运行结果：
the math sol is [0.09090909 0.63636364]
[array([-8, -3])]
[0.22054380664652568]
[array([-8, -3])]
the numerical sol is [array([2, 1]), array([0.23564955, 0.33836858]), array([0.09090909, 0.63636364])]


展开全文
• 拟牛顿法、高斯牛顿法、牛顿法、共轭梯度法法的python实现，使用的测试案例都来自于课本例题便于验证。

本文的代码与课本算法的流程框架完全一致，不再重复说明理论，只是汇总一下代码实现。

文章目录
拟牛顿法-BFGS高斯-GN牛顿法-Newton共轭梯度法总结

拟牛顿法-BFGS
# -*- coding: utf-8 -*-#
# Author: xhc
# Date:    2021-05-28 16:01
# project:  1_zyh
# Name:    BFGS.py
import numpy as np

# 定义求解方程
fun = lambda x: 0.5*x[0]**2+x[1]**2-x[0]*x[1]-2*x[0]

# 定义函数 用于求梯度数值 数据放入mat矩阵
def gfun(x):
x = np.array(x)
part1 = x[0][0]-x[1][0]-2
part2 = -1.0*x[0][0]+2*x[1][0]
result = np.mat([part1,part2])
# print(gf) 打印测试
return result

# 传入参数有三个 fun为函数 gfun为梯度 x0为初始值
def BGFS(fun, gfun, x0):
maxk = 5000 # 迭代次数
rho = 0.45 # 步长
sigma = 0.3 # 常数 在 0到1/2之间
epsilon = 1e-6 # 终止条件
k = 0 # 第几次迭代
Bk = np.mat([[1.0,0.0],[0.0,1.0]]) # 初始值为单位矩阵

while k < maxk:
gk = gfun(x0) # (1,2)
if np.linalg.norm(gk) < epsilon: # 范数小于epsilon 终止
break
dk = -Bk.I.dot(gk.T) # Bk*d+梯度=0
m = 0
mk = 0
while m < 30: # Arijo算法求步长
if fun(x0 + (rho**m)*dk) < fun(x0) + sigma * (rho**m)*gk.dot(dk):
mk = m
break
m = m + 1

x = x0 + (rho**mk)*dk # (2,1) # 迭代公式
# 以下就是BFGS修正公式
sk = x - x0 # (2,1)
yk = gfun(x) - gk  # (1,2)
if yk * sk > 0:
Bk = Bk - (Bk * sk * sk.T * Bk)/(sk.T*Bk*sk) + (yk.T*yk)/(yk*sk)
k = k + 1
x0 = x
print("--------------------------")
print("当前点为为%s" % x0.T)
print("当前点的值为%f" % fun(x0))
print("--------------------------")

return x, k
if __name__ == "__main__":
x0 = np.array([[1],[1]]) # 选择初始点
x,  k = BGFS(fun, gfun, x0)
print("BFGS")
print("本题使用的例子是课本p56例4.1.1，但是采用的非精确搜索，结果如下")
print("迭代次数为%d次" %k)
print("最优点为%s" %x.T)
print("最小值为%d" %fun(x))



高斯-GN
# -*- coding: utf-8 -*-#
# Author: xhc
# Date:    2021-05-29 14:56
# project:  1_zyh
# Name:    GN.py

import numpy as np
#fun =

def gfun(x):
x = np.array(x)
val1 = x[0][0] - 0.7*np.sin(x[0][0]) - 0.2*np.cos(x[1][0])
val2 = x[1][0] - 0.7*np.cos(x[0][0]) + 0.2*np.sin(x[1][0])
y = np.mat([[val1],[val2]])
return y

def Hess(x):
x = np.array(x)
val1 = 1 - 0.7*np.cos(x[0][0])
val2 = 0.2*np.sin(x[1][0])
val3 = 0.7*np.sin(x[0][0])
val4 = 1 + 0.2*np.cos(x[1][0])
y=np.mat([[val1, val2],
[val3, val4]])
return y

def GN(gfun, Hess, x0):
maxk = 5000
rho = 0.4
sigma = 0.4
k = 0
epsilon = 1e-6
while k < maxk:
fk = gfun(x0)
hess = Hess(x0)

gk = hess.T * fk

dk = - (hess.T * hess ).I.dot(gk)

if np.linalg.norm(gk) < epsilon:
break
m = 0
mk = 0
while m < 30:
newf = 0.5 * np.linalg.norm(gfun(x0 + (rho**m) * dk))**2
oldf = 0.5 * np.linalg.norm(gfun(x0))**2
if newf < oldf + sigma * (rho**m)*gk.T*dk:
mk = m
break
m = m + 1

x0 = x0 + (rho**mk) * dk
k = k + 1
print("--------------------------")
print("当前点为为%s" % x0.T)
print("--------------------------")
x = x0
return x,  k
if __name__ == "__main__":
x0 = np.array([[0], [0]])  # (2,1)
x, k = GN(gfun, Hess, x0)
print("GN")
print("本题实现非独立完成，在课本找不到例题，最初使用1/2(cos^2(x1)+1/2(sin^2(x2)) 会产生奇异矩阵，无法求逆。更换函数 结果如下")
print("迭代次数为%d次" % k)
print("最优点为%s" % x.T)
# print("最小值为%d" %fun(x))



牛顿法-Newton
# -*- coding: utf-8 -*-#
# Author: xhc
# Date:    2021-05-28 17:06
# project:  1_zyh
# Name:    Newton.py
import numpy as np

fun = lambda x:0.5*x[0]**2+x[1]**2-x[0]*x[1]-x[0]

def gfun(x):
x = np.array(x)
part1 = x[0][0] - x[1][0] - 1
part2 = -x[0][0] + 2 * x[1][0]
gf = np.mat([part1,part2])
return gf

def Hess(x):
x = np.array(x)
part1 = 1
part2 = -1
part3 = -1
part4 = 2
Hess = np.mat([[part1, part2],[part3, part4]])

return Hess

def DampNewtonMethod(fun, gfun, Hessian, x0):
maxk = 5000 # 最大迭代次数
rho = 0.4 # 步长
sigma = 0.4 # 常数
k = 0 # 迭代次数
epsilon = 1e-6 # 终止条件
while k < maxk:
gk = gfun(x0) # 梯度计算
Gk = Hessian(x0) # hess计算
dk = -Gk.I.dot(gk.T) # 搜索方向
if np.linalg.norm(gk) < epsilon: # 终止条件
break
m = 0
mk = 0
while m < 20: # 用Armijo搜索求步长
if fun(x0 + (rho**m) * dk) < fun(x0) + (sigma * (rho**m) * gk.dot(dk)):
mk = m
break
m = m + 1
x0 = x0 + (rho**mk) * dk
k = k + 1
print("--------------------------")
print("当前点为为%s" % x0.T)
print("当前点的值为%f" % fun(x0))
print("--------------------------")
x = x0
return x, k
if __name__ == "__main__":
x0 = np.array([[0],[0]])
x, k = DampNewtonMethod(fun, gfun, Hess, x0)
print("牛顿法")
print("本题使用的例子是课本p41例3.2.1，但是采用的非精确搜索，结果如下")
print("迭代次数为%d次" %k)
print("最优点为%s" %x.T)
print("最小值为%f" %fun(x))



共轭梯度法
# -*- coding: utf-8 -*-#
# Author: xhc
# Date:    2021-05-29 15:40
# project:  1_zyh
# Name:    共轭梯度法法.py

import numpy as np

# 函数
fun = lambda x: 0.5*x[0]**2+x[1]**2

# 梯度
def gfun(x):
x = np.array(x)
part1 = x[0][0]
part2 = 2*x[1][0]
gf = np.mat([part1,part2])
return gf

maxk = 5000
rho = 0.4
sigma = 0.4
k = 0
epsilon = 1e-6
n = len(x0)
g0 = 0
d0 = 0
while k < maxk:
g = gfun(x0)
itern = k - (n+1)*np.floor(k/(n+1))
itern = itern + 1
if itern == 1: # 每一次迭代k=0的情况 负梯度
d = -g
else: # 每一次迭代 k>0 的情况  负梯度+bkd
beta = (g.dot(g.T))/(g0.dot(g0.T))  # 更新beta
d = -g + beta*d0
gd = g.dot(d.T)
if gd >= 0.0:
d = -g
if np.linalg.norm(g) < epsilon: # 终止条件
break
m = 0
mk = 0
# armijo
while m < 20:
if fun(x0 + (rho**m)*d.T) < fun(x0) + sigma * (rho**m) * g.dot(d.T):
mk = m
break
m = m + 1
x0 = x0 + (rho**mk)*d.T
x0 = np.array(x0)
g0 = g
d0 = d
k = k + 1
print("--------------------------")
print("当前点为为%s" % x0.T)
print("当前点的值为%f" % fun(x0))
print("--------------------------")
x = x0
return x, k
if __name__ == "__main__":
x0 = np.array([[2],[1]])
x, k = FR_Gradient(fun, gfun, x0)
print("共轭梯度法--FR算法")
print("本题使用的例子是课本p76例5.2.1")
print("迭代次数为%d次" %k)
print("最优点为%s" %x.T)
print("最小值为%f"%fun(x))



总结

如有错误，欢迎指出。


展开全文
• 共轭梯度法也是解决无约束优化问题的常用迭代算法，它结合了最速下降法矩阵共轭梯度的性质，可以加快算法的迭代过程。且如果初始点选取后的最终优化中不满足精度条件，还可保存上一步得到的迭代点进行再次迭代直到...
• 共轭梯度法（以希尔伯特矩阵为例） 实现共轭梯度算法并且使用它解决对称矩阵——希尔伯特矩阵为系数矩阵A时，右边常数项列向量为单位列向量，初始点为零向量。展示维度分别为5，8，12，20时的迭代次数与结果。 迭代...
• main为主函数 fun gfun ggfun分别为输入的...分别为最速下降法 牛顿法（阻尼）共轭梯度法 以及 拟牛顿法 F1-4为下降的图示 可以看到牛顿法和拟牛顿法收敛速度最快 但是牛顿法需要求矩阵的逆 在实际中 运算量可能较大
• python实现共轭梯度
• 共轭梯度法求解线性方程组 Ax = b 的子程序，其中，A为系数矩阵，x为待求未知数向量，b为方程右端向量。 具体算法见 https://wenku.baidu.com/view/f5157b0c9a6648d7c1c708a1284ac850ac02047c.html?rec_flag=default...
• 待写
• % 计算xk点的梯度梯度值 fun=fun(x1,x2); fx1=diff(fun,'x1'); fx2=diff(fun,'x2'); fun=inline(fun); fx1=inline(fx1); fx2=inline(fx2); funval=feval(fun,xk1(1),xk1(2)); gradx1=feval(fx1,xk1(1&...
• 今天小编就为大家分享一篇python 梯度法求解函数极值的实例，具有很好的参考价值，希望对大家有所帮助。一起跟随小编过来看看吧
• 西安交通大学《数值分析》第三章课后题3.2 import numpy as np import math from math import log import matplotlib.pyplot as plt def generate_... 1.] count 25 另一篇笔记 快速下降需要计算快8k次，我人没了。
• 程序员的数学-优化理论
• Python实现共轭梯度法与最速下降法的比较
• 共轭梯度法也是共轭方向法中的一种，但是它减少了梯度方向的搜索量，它直接采取经过一维搜索最小点处的梯度方向作为我们的搜索方向，因而在计算速度上有了一定的提升。如果你对这些优化算法感到困惑，现在你需要明白...
• 个人博客：blog.moon.top
• 梯度法求解线性方程组：Ax=b 数值分析老师给的作业，写出来代码平时分满分。简单勿喷，hhhh 原理 CG求解线性方程组的原理大家翻一翻数值分析的课本即可，此处不哔哔直接上matlab代码 代码 %qjyang clc clear; A = ...
• ## 共轭梯度法原理与实现

万次阅读 多人点赞 2015-05-29 23:24:47
作者：金良（golden1314521@gmail.com） csdn博客： 共轭...定义 共轭方向的性质 共轭方向法 算法描述 算法的收敛性 搜索步长kalpha_k的确定 共轭梯度法 共轭梯度法的原理 共轭梯度算法描述 共轭梯度算法Python实现 r

...

python 订阅