Repository navigation
Expand file tree
/
Copy pathBFGS.py
More file actions
75 lines (57 loc) · 2 KB
/
Copy pathBFGS.py
File metadata and controls
75 lines (57 loc) · 2 KB
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
from mcholmz import modifiedChol
import numpy as np
def inexact_line_search(x, dir, func, grad):
sigma = 0.25
beta = 0.5
a0 = 1
f0 = func(x)
# s, phi1 = ls.scalar_search_armijo(func, x, d, f0, alpha0=1)
c = (grad(x).T @ dir) # =phi(0)
return wolf(x, dir, sigma, 1, a0, func, grad, f0, beta, c)
def wolf(x, d, sigma, b, a0, f, g, f0, cb, c):
alpha = a0 * b
bound = (sigma * alpha * c) + f0
# bound = (g(x) @ d)*a*b
x_ad = x + alpha * d
phi = f(x_ad)
g2 = g(x_ad).T @ d
g1 = 0.9 * c
if bound >= phi and g2 >= g1: # 0.9 is c2 this is wolfe curvature condition
# if bound >= phi and g(x + alpha * d) @ d >= 0.9 * c: # 0.9 is c2 this is wolfe curvature condition
# if s * bound >= phi:
return alpha
return wolf(x, d, sigma, b * cb, a0, f, g, f0, cb, c)
def update(Bk, pk, qk):
pk = pk.reshape(pk.shape[0], 1)
qk = qk.reshape(pk.shape[0], 1)
# Broyden family
sk = np.matmul(Bk, qk)
tk = sk.T @ qk
miuk = pk.T @ qk
vk = (1 / miuk) * pk - (1 / tk) * sk
norm = np.linalg.norm
return Bk + (1 / miuk) * pk @ pk.T - (1 / tk) * sk @ sk.T + tk * vk @ vk.T # return Bk+1
def BFGS(func, gradient, start_point, learn_rate, exact_ls=True, stop_criteria=10 ** -5):
zero = 10 ** -10
x = np.array(start_point)
i = 0
B = np.eye(x.shape[0]) # B0 initialization
g = -gradient(x)
tod = np.array([start_point])
prev_x = x
prev_g = g
while np.linalg.norm(g) > stop_criteria:
i += 1
dk = B @ g
if exact_ls is True:
a = learn_rate(x, dk)
else:
a = learn_rate(x, dk, func, gradient)
prev_x = np.array(x) # keep old x and g for p and q accordingly
prev_g = np.array(g)
x += dk * a # calc new x and g
tod = np.append(tod, [x], axis=0)
g = -gradient(x)
# print(f"{i}:{x}")
B = update(B, x - prev_x, -(g - prev_g)) # update approx inv hessian
return i, tod # returns optimal x