I want to implement a conjugate gradient algorithm in python. However, when I run the code no results are printed/displayed. When I stop the program python shows me this:
def V(r):
return (r[0]-r[1])**4+2*r[0]**2+r[1]**2-r[0]+2*r[1]
def g(r):
return np.array([4*(r[0]-r[1])**3+4*r[0]-1, -4*(r[0]-r[1])**3+2*r[1]+2])
def Conj_Grad(r):
iterations = 0
maxiterations = 1000
precision = 10**-10
a = 10**-4
b = 0.9
alpha = 0.09
r0 = r
g0 = g(r)
p0 = -g(r)
V0 = V(r)
Coord_List= []
S = np.dot(p0,g0)
while iterations < maxiterations:
Coord_List.append(r0)
while True:
r1 = r0 + alpha*p0
g1 = g(r1)
V1 = V(r1)
if V1 <= V0 + a*alpha*S:
g1 = g(r1)
if abs(np.dot(g1,p0)) <= b*abs(np.dot(p0,g0)):
break
else:
alpha = alpha*(1.5)
continue
else:
alpha = alpha/2
continue
if abs(V0-V1) < precision:
Coord_List.append(r1)
iterations +=1
break
# Generating beta by Fletcher Reeves
beta = np.dot(g1,g1)/np.dot(g0,g0)
p0 = -g1 + beta*p0
g0 = g1
V0 = V1
iterations +=1
return (Coord_List, i, V1-V0)
<ipython-input-57-f26ac2ab8296> in Conj_Grad(r)
29 Coord_List.append(r0)
30 while True:
---> 31 r1 = r0 + alpha*p0
32 g1 = g(r1)
33 V1 = V(r1)
When I input an r, for example [1,1]; no matter how long or short I let it run it always seems to be stuck at line 31. What is wrong with this code?
while True
is an infinite loop. See how to create a minimal reproducible example.