added backtracking to: TIRUNEH, Ababu Teklemariam. A modified three-point Secant method with improved rate and characteristics of convergence. arXiv preprint arXiv:1902.09058, 2019.

3psecant.py

import matplotlib.pyplot as plt
from matplotlib import patches
from matplotlib.collections import LineCollection
from numpy import finfo, array, sqrt, isnan, concatenate, sin, exp


def secant_ls_3p(y, x_0, tol, max_it=100):
    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
    rellength = abs(p / x_k)
    lambda_min = tol / rellength
    accum_step = 0.0
    total_backtracks = 0
    backtrackcount = 0
    steps = []
    y_list = []
    x_list = []
    for j in range(max_it):
        y_list += [y_k]
        x_list += [x_k]
        steps += [accum_step]
        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))
        elif 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.4g}' +
                   '\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:
            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)
                    )
                )
            print(('{:d}:\t x_k: {:0.4f}\ty_k: {:0.4g}' +
                   '\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

        stop = False
        lambda_ls = 1.0
        total_backtracks += backtrackcount
        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
            accum_step += lambda_ls
            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

            x_list += [x_2]
            y_list += [f_2]
            steps += [accum_step]

            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
                backtrackcount += 1
                accum_step -= lambda_ls
                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
        if abs(p) <= tol:
            # avoid 1/0 division
            break

    print(('{:d}:\t x_k: {:0.4f}\ty_k: {:0.4g}'+
                    '\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))

    soln = dict()
    for item in ['x_k', 'y_k', 'x_list', 'y_list',
                 'j', 'total_backtracks', 'steps']:
        soln[item] = locals().get(item)
    return soln


def plot_results(x_list, y_list, steps, j, total_backtracks, fun_title=''):
    fig, axes = plt.subplots(1, 2, sharey=False)
    axes[0].plot(x_list, y_list, 'o', fillstyle='none', markeredgewidth=0.5)
    for i in range(len(x_list) - 1):
        if not isnan(x_list[i]) and not isnan(x_list[i + 1]):
            dx = x_list[i + 1] - x_list[i]
            dy = y_list[i + 1] - y_list[i]
            arrow_tail = [x_list[i], y_list[i]]
            arrow_head = [x_list[i + 1], y_list[i + 1]]
            axes[0].arrow(
                x_list[i], y_list[i], dx, dy, head_width=0.02, head_length=0.2, fc='k', ec='k'
            )

    x_list_not_nan = [x for x in x_list if not isnan(x)]
    y_list_not_nan = [x for x in y_list if not isnan(x)]
    i_list_not_nan = range(len(x_list_not_nan))
    points = array([x_list_not_nan, y_list_not_nan]).T.reshape(-1, 1, 2)
    segments = concatenate([points[:-1], points[1:]], axis=1)
    norm = plt.Normalize(0, len(x_list_not_nan))
    lc = LineCollection(segments, cmap='viridis', norm=norm)

    ax, fig = axes[0], plt.gcf()
    lc.set_array(array(i_list_not_nan))
    lc.set_linewidth(2)
    line = ax.add_collection(lc)
    axes[0].set_xlabel('x')
    axes[0].set_ylabel('f')
    axes[0].set_title('$'+fun_title.replace('**', '^')+'$')

    axes[1].loglog(array(steps), abs(array(y_list)))
    axes[1].set_xlabel('steps - accum.')
    axes[1].set_ylabel('abs(f)')
    axes[1].set_title('iterations: {:d}\n tot. backtracks: {:d}'.format(
        j, total_backtracks))

    cb = fig.colorbar(line)
    cb.set_label('iterations')
    plt.tight_layout()

tol = finfo(float).eps
alpha = 1e-4


f1 = 'x**3 + 4 * x**2 - 10'
x_0, x_1, x_2 = 0.5, 0.55, 0.6
soln = secant_ls_3p(lambda x: eval(f1), x_0, tol)
plot_results(soln['x_list'], soln['y_list'], soln['steps'], soln['j'], soln['total_backtracks'], fun_title=f1)
print(5*'\n')

x_0, x_1, x_2 = 0.9, 0.95, 1.0
soln = secant_ls_3p(lambda x: eval(f1), x_0, tol)
plot_results(soln['x_list'], soln['y_list'], soln['steps'], soln['j'], soln['total_backtracks'], fun_title=f1)
print(5*'\n')

f2 = 'sin(x)**2 - x**2 + 1'
x_0 = -1.0
soln = secant_ls_3p(lambda x: eval(f2), x_0, tol)
plot_results(soln['x_list'], soln['y_list'], soln['steps'], soln['j'], soln['total_backtracks'], fun_title=f2)
print(5*'\n')

x_0 = -3.5
soln = secant_ls_3p(lambda x: eval(f2), x_0, tol)
plot_results(soln['x_list'], soln['y_list'], soln['steps'], soln['j'], soln['total_backtracks'], fun_title=f2)
print(5*'\n')

f3 = '(x - 2) * (x + 2)**4'
x_0 = -3.1
soln = secant_ls_3p(lambda x: eval(f3), x_0, tol)
plot_results(soln['x_list'], soln['y_list'], soln['steps'], soln['j'], soln['total_backtracks'], fun_title=f3)
print(5*'\n')

f4 = 'exp(x**2+7*x-30)-1'
x_0 = 4.0
soln = secant_ls_3p(lambda x: eval(f4), x_0, tol)
plot_results(soln['x_list'], soln['y_list'], soln['steps'], soln['j'], soln['total_backtracks'], fun_title=f4)
print(5*'\n')


f5 = '(x - 1)**6 - 1'
x_0 = 1.5
soln = secant_ls_3p(lambda x: eval(f5), x_0, tol)
plot_results(soln['x_list'], soln['y_list'], soln['steps'], soln['j'], soln['total_backtracks'], fun_title=f5)
print(5*'\n')

plt.show()