How can I deploy a CG-Steihaug algorithm for trust region sub-problem solving?

Question

I'm studying various optimization methods and on this occasion, I'm trying to tackle the Trust Region Problem by solving the sub-region problem with the Steihaug-CG algorithm in Python. I'm using the pseudo-code from the book (Numerical Optimization, Nocedal), which is the same Pseudo Code from this source. I'm more than confused on the both find $\tau$ steps.

 

I tried several ways to decompose this, mainly by using a solver for the $m_k{(p_k)}$ function and a line search for the second possible case, but both proved completely unsuccessful. I'm also adding my Python code for further clarification:

def steihaugcg(B, gradf, delta, tol=1e-9, max_it=1000):
r=[gradf]
if norm(r[-1]) < tol: return np.zeros(B.shape[0])

def LineSearch(z, d, DELTA):
    t=.5
    for _ in range(500):
        if np.allclose(norm(z+t*d),DELTA):return t
        if norm(z+t*d) < DELTA: t = t*1.9
        if norm(z+t*d) > DELTA: t = t*0.1
    return None

d=-r[-1]
t=0
Size=B.shape[0]
z=np.zeros(Size)
for _ in range(max_it):
    if d.T@B@d <= 0:
        t = minimize(lambda t: gradf.T@(z+t*d) + 0.5*(z+t*d).T@B@(z+t*d), 1).x
        p=z+t*d
        if np.allclose(norm(p),delta):
            return p
    alpha=(r[-1].T@r[-1])/(d.T@B@d)
    z=z+alpha*d

    if norm(z) >= delta:
        t = LineSearch(z,d,delta)
        if t is not None: return z+t*d

    r.append(r[-1]+alpha*B@d)

    if norm(r[-1])<tol:print("third");return z
    beta = (r[-1].T@r[-1])/(r[-2].T@r[-2])
    d = -r[-1] + beta*d

Can somebody provide insight on how to solve the find $\tau$ subproblem?

Answer 1

Note that $\|p_k\|^2=\|z_j+\tau d_j\|^2 = z_j^Tz_j+2\tau d_j^Tz_j+\tau^2d_j^Td_j=\Delta_k^2$. So you can solve $\tau$ using the quadratic formula.