backtracking secant method based on: Jr., J. E. Dennis ; Schnabel, Robert B.: Numerical Methods for Unconstrained Optimization and Nonlinear Equations. Philadelphia: SIAM, 1996.

secant_backtrack.py

import matplotlib.pyplot as plt
from numpy import finfo, array, sqrt

def y(x):
    return x**3 + 4 * x**2 - 10

max_it = 100
tol = finfo(float).eps
alpha = 1e-4
y_list = []
x_list = []
x_0, x_1, x_2 = 0.5, 0.55, 0.6
# x_0, x_1, x_2 = 0.9, 0.95, 1.0
# p = x_1 - x_0
x_k = x_0
y_k = y(x_k)
g_k = 1 / 2 * y_k**2
g_prime_k = - y_k**2
# relative length of p as calculated in the stopping routine
# p = -0.1 * x_0
p = (x_1 - x_0)
rellength = abs(p / x_k)
lambda_min = tol / rellength
for j in range(max_it):
    y_list += [y_k]
    x_list += [x_k]
    if abs(y_k) <= tol:
        break
    if j == 0:
        inv_slope = 0.1 / y_k
        p = -inv_slope * y_k
        x_k_minus_1 = x_k
        y_k_minus_1 = y_k
        g_k_minus_1 = g_k
        g_prime_k_minus_1 = g_prime_k

        print(('{:d}:\t x_k: {:0.4f}\ty_k: {:0.4f}' +
            '\tx_k_minus_1: {:0.4f}'
            ).format(j, x_k, y_k, x_k_minus_1))
    else:#if j == 1:
        inv_slope = (x_k - x_k_minus_1) / (y_k - y_k_minus_1)
        # x_k_plus_1 = x_k - inv_slope * y_k
        p = -inv_slope * y_k
        # x_k_plus_1 = x_2
        # p = (x_2 - x_1)
        print(('{:d}:\t x_k: {:0.4f}\ty_k: {:0.4f}' +
               '\tx_k_minus_1: {:0.4f}, p: {:0.4f}'
               ).format(j, x_k, y_k, x_k_minus_1, p))

        x_k_minus_2 = x_k_minus_1
        x_k_minus_1 = x_k

        y_k_minus_2 = y_k_minus_1
        y_k_minus_1 = y_k
        g_k_minus_1 = g_k
        g_prime_k_minus_1 = g_prime_k
    # else:
    #     # x_k_plus_1 = x_k_minus_2 - y_k_minus_2 * (
    #     #     y_k - y_k_minus_1
    #     #     ) / ((y_k - y_k_minus_2) / (x_k - x_k_minus_2) * (
    #     #         y_k - y_k_minus_1) - y_k * (
    #     #         (y_k - y_k_minus_2) / (x_k - x_k_minus_2) -
    #     #         (y_k_minus_1 - y_k_minus_2) / (x_k_minus_1 - x_k_minus_2)
    #     #         )
    #     #     )
    #     p = (x_k_minus_2 - x_k) - y_k_minus_2 * (
    #         y_k - y_k_minus_1
    #         ) / ((y_k - y_k_minus_2) / (x_k - x_k_minus_2) * (
    #             y_k - y_k_minus_1) - y_k * (
    #             (y_k - y_k_minus_2) / (x_k - x_k_minus_2) -
    #             (y_k_minus_1 - y_k_minus_2) / (x_k_minus_1 - x_k_minus_2)
    #             )
    #         )
    #     x_k_minus_2 = x_k_minus_1
    #     x_k_minus_1 = x_k
    #
    #     y_k_minus_2 = y_k_minus_1
    #     y_k_minus_1 = y_k
    #     g_k_minus_1 = g_k
    #     g_prime_k_minus_1 = g_prime_k
    #     print(('{:d}:\t x_k: {:0.4f}\ty_k: {:0.4f}'+
    #         '\tx_k_minus_1: {:0.4f}\tx_k_minus_2: {:0.4f}'
    #         ).format(j, x_k, y_k, x_k_minus_1, x_k_minus_2))
    
    stop = False
    lambda_ls = 1.0
    backtrackcount = 0
    g_0 = g_k_minus_1
    g_prime_0 = g_prime_k_minus_1
    f_0 = y_k_minus_1
    while not stop and j + backtrackcount <= max_it:
        # backtracking, line search - numerical recipes 3ed
        backtrackcount += 1
        if lambda_ls <= lambda_min:
            # opposite direction
            p = -p
            lambda_ls = 1.0
            stop = False

        x_2 = x_k + lambda_ls * p
        f_2 = y(x_2)
        g_2 = 1 / 2 * f_2**2
        descent = alpha * lambda_ls * g_prime_0
        g_max = g_0 + descent

        satisfactory = g_2 <= g_max
        stop = satisfactory or lambda_ls <= lambda_min
        print(('{:d}-{:d}:\t x_2: {:.2g}\ty_2: {:.2g}' +
               '\tg_0: {:.2g}\tg_2: {:.2g}\tg_max: {:.2g}\tlambda: {:.2g}\t p: {:2g}').format(
            j, backtrackcount, x_2, f_2, g_0, g_2, g_max, lambda_ls, p))

        if not stop:
            # backtrack - reduce lambda
            if lambda_ls == 1:
                # first backtrack quadratic fit
                lambda_temp = -g_prime_0 / (
                        2 * (g_2 - g_0 - g_prime_0)
                    )
            elif lambda_ls < 1:
                # subsequent backtracks cubic fit
                a, b = 1 / (lambda_ls - lambda_prev) * array(
                        [[+1 / lambda_ls**2, -1 / lambda_prev**2],
                        [-lambda_prev / lambda_ls**2, 
                         +lambda_ls / lambda_prev**2]]
                    ).dot(array(
                        [[g_2 - g_0 - g_prime_0 * lambda_ls],
                        [ g_1 - g_0 - g_prime_0 * lambda_prev]]
                        ))
                a, b = a.item(), b.item()
                disc = b**2 - 3 * a * g_prime_0
                if a == 0:
                    # actually quadratic
                    lambda_temp = - g_prime_0 / (2 * b)
                else:
                    # legitimate cubic
                    lambda_temp = (-b + sqrt(disc)) / (3 * a)
                if lambda_temp > 1 / 2 * lambda_ls:
                    lambda_temp = 1 / 2 * lambda_ls                    
            lambda_prev = lambda_ls
            g_1 = g_2
            if lambda_temp <= 0.1 * lambda_ls:
                lambda_ls = 0.1 * lambda_ls
            else:
                lambda_ls = lambda_temp
        x_k_plus_1 = x_2
        y_k_plus_1 = f_2
        g_k_plus_1 = g_2


    x_k = x_k_plus_1
    y_k = y_k_plus_1
    g_k = g_k_plus_1
    g_prime_k = - y_k ** 2

print(('{:d}:\t x_k: {:0.4f}\ty_k: {:0.4f}'+
                '\tx_k_minus_1: {:0.4f}\tx_k_minus_2: {:0.4f}'
                ).format(j, x_k, y_k, x_k_minus_1, x_k_minus_2))