服务器之家

服务器之家 > 正文

python实现共轭梯度法

时间:2021-07-30 00:53     来源/作者:Tomator01

共轭梯度法是介于最速下降法与牛顿法之间的一个方法,它仅需利用一阶导数信息,但克服了最速下降法收敛慢的缺点,又避免了牛顿法需要存储和计算Hesse矩阵并求逆的缺点,共轭梯度法不仅是解决大型线性方程组最有用的方法之一,也是解大型非线性最优化最有效的算法之一。 在各种优化算法中,共轭梯度法是非常重要的一种。其优点是所需存储量小,具有步收敛性,稳定性高,而且不需要任何外来参数。

算法步骤:

python实现共轭梯度法

?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
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是试验的部分,在本段程序中不起作用。

迭代轨迹:

python实现共轭梯度法

python实现共轭梯度法

三种最优化方法的迭代次数对比:

 

最优化方法

最速下降法

共轭梯度法

牛顿法

迭代次数

1702

240

5

 

以上就是本文的全部内容,希望对大家的学习有所帮助,也希望大家多多支持服务器之家。

原文链接:https://blog.csdn.net/Big_Pai/article/details/88540147

相关文章

热门资讯

2020微信伤感网名听哭了 让对方看到心疼的伤感网名大全
2020微信伤感网名听哭了 让对方看到心疼的伤感网名大全 2019-12-26
yue是什么意思 网络流行语yue了是什么梗
yue是什么意思 网络流行语yue了是什么梗 2020-10-11
背刺什么意思 网络词语背刺是什么梗
背刺什么意思 网络词语背刺是什么梗 2020-05-22
苹果12mini价格表官网报价 iPhone12mini全版本价格汇总
苹果12mini价格表官网报价 iPhone12mini全版本价格汇总 2020-11-13
2021德云社封箱演出完整版 2021年德云社封箱演出在线看
2021德云社封箱演出完整版 2021年德云社封箱演出在线看 2021-03-15
返回顶部