Goldstein準則優化演算法（julia實現）

來自專欄演算法julia實現練習

"""

計算向量範數

1）給定向量x，範數p（預設2）

2）(Σi xi^p)^(1/p)

"""

function norm(x,p=2)

sum(x.^p)^(1/p)

end

"""

計算搜索步長

1）給定當前點x，當前函數值y，方向d，當前步長a，梯度長度gl，初始化搜索範圍(a0,b0)=(0,inf)

2）計算新點x1=x+ad，函數值y1=f(x1)

3）如果新點滿足Armijo、Goldstein條件，演算法終止

4）如果新點不滿足Armijo條件，修改搜索範圍上限b0=a，縮小a=a0+(a-a0)*s

5）如果新點不滿足Goldstein條件，修改搜索範圍下限a0=a，放大a=a/s

6）轉2

"""

function findstep(x, y, f, d, a, gl; r=0.001, s=0.618, maxiter=10000)

a0 = 0

b0 = Inf

x1 = x + a*d

y1 = f(x1)

for k = 1:maxiter

if y1 > y - r*gl*a

b0 = a

a = a0 + (a-a0)*s

elseif y1 < y - (1-r)*gl*a

b0 = a

a = a/s

else

break

end

x1 = x + a*d

y1 = f(x1)

end

a, x1, y1

end

"""

Goldstein準則的無約束最優化演算法

1）給定初始點x=x0，初始步長a=a0

2）計算函數值y=f(x)，負梯度方向d=-g(x)

3）如果負梯度長度很小，演算法終止

4）尋找合適步長，滿足Goldstein準則，計算新步長a、更新點x、函數值y

5）判斷函數值變化是否很小，如果很小演算法終止，如果達到迭代最大次數演算法終止，否則轉2

"""

function optim(x0, f, g; maxiter=10000, e=1e-6, a0=0.001)

x = x0

a = a0

for k = 1:maxiter

y = f(x)

grad = g(x)

gl = norm(grad)

if gl^0.5 < e

break

end

d = -grad / gl

a, x1, y1 = findstep(x, y, f, d, a, gl)

if k % 10 == 0

println("k=", k, ", y1=", y1, ", x=", x, ", a=", a)

end

if maximum(abs.(x1 .- x)) < e

break

end

x = x1

end

x

end

#測試

function test()

function f(x)

sum(x.^2)

end

function g(x)

2*x

end

optim([1,2,3], f, g)

end

#reference

# 《數字最優化方法》高立 演算法2.1 Armijo準則 Goldstein準則

推薦閱讀：