python實現共軛梯度法

2022-10-04 17:12:14 字數 3840 閱讀 2753

共軛梯度法是介於最速下降法與牛頓法之間的乙個方法,它僅需利用一階導數資訊,但克服了最速下降法收斂慢的缺點,又避免了牛頓法需要儲存和計算hesse矩陣並求逆的缺點,共軛梯度法不僅是解決大型線性方程組最www.cppcns.com有用的方法之一,也是解大型非線性最優化最有效的演算法之一。 在各種優化演算法中,共軛梯度法是非常重要的一種。其優點是所需儲存量小,具有步收斂性,穩定性高,而且不需要任何外來引數。

演算法步驟:

import random

import numpy as np

import matplotlib.pyplot as plt

def goldsteinsearch(f,df,d,x,alpham,rho,t):

'''線性搜尋子函式

數f,導數df,當前迭代點x和當前搜尋方向d,t1,

'''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是試驗的部分,在本段程式中不起作用。

迭代軌跡:

三種最優化方法的迭代次數對比:

最優化方法

最速下降法

共軛梯度法

牛頓法程式設計客棧

迭代次數

www.cppcns.com

程式設計客棧

本文標題: python實現共軛梯度法

本文位址:

python實現共軛梯度法

共軛梯度法 conjugate gradient 是介於最速下降法與牛頓法之間的乙個方法,它僅需利用一階導數資訊,但克服了最速下降法收斂慢的缺點,又避免了牛頓法需要儲存和計算hesse矩陣並求逆的缺點,共軛梯度法不僅是解決大型線性方程組最有用的方法之一,也是解大型非線性最優化最有效的演算法之一。這裡...

21 梯度法及共軛梯度法

coding utf 8 created on wed mar 20 21 56 21 2019 author zhangchaoyu 當離最優值較遠時,利用梯度法,接近極小點時用收斂較快的其他方法 import math import copy 1 f x,y x x 25 y y 2 f x,y...

用MATLAB實現FR共軛梯度法求解例項

問題 編寫fr共軛梯度法用於求解min x2 1 x1 x2 x 22 2 x1 4 x2初始點取為 x0 2,2 t frcg.m function x,val,k frcg fun,funs,x0 功能 用fr共軛梯度法求解無約束問題 min f x 輸入 x0是初始點,fun,gfun分別是目...