optimize.py

#Trust Region (Dogleg) vs Steepest Descent

import numpy as np
import matplotlib.pyplot as plt

# Problem Setup
Q = np.array([
    [100, 20, 0, 0, 0],
    [20, 50, 10, 0, 0],
    [0, 10, 20, 5, 0],
    [0, 0, 5, 10, 2],
    [0, 0, 0, 2, 5]
], dtype=float)

b = np.array([1, 2, 3, 4, 5], dtype=float)
x0 = np.zeros(5)

def f(x):
    return 0.5 * x.T @ Q @ x - b.T @ x

def grad_f(x):
    return Q @ x - b

# Steepest Descent with Exact Line Search
def steepest_descent(x0, tol=1e-8, maxiter=200):
    x = x0.copy()
    hist = {'x': [], 'f': [], 'g': []}
    
    for k in range(maxiter):
        g = grad_f(x)
        hist['x'].append(x.copy())
        hist['f'].append(f(x))
        hist['g'].append(np.linalg.norm(g))
        
        if np.linalg.norm(g) < tol:
            break
        
        alpha = (g.T @ g) / (g.T @ Q @ g)
        x = x - alpha * g
    
    return hist, x

# Dogleg Trust Region Method
def trust_region_dogleg(x0, Delta0=1.0, Delta_max=100.0, eta=0.15, tol=1e-8, maxiter=200):
    x = x0.copy()
    Delta = Delta0
    hist = {'x': [], 'f': [], 'g': [], 'Delta': []}
    
    for k in range(maxiter):
        g = grad_f(x)
        hist['x'].append(x.copy())
        hist['f'].append(f(x))
        hist['g'].append(np.linalg.norm(g))
        hist['Delta'].append(Delta)
        
        if np.linalg.norm(g) < tol:
            break
        
        # Newton and Cauchy steps
        pB = -np.linalg.solve(Q, g)
        gBg = g.T @ Q @ g
        pU = -(g.T @ g) / gBg * g
        
        if np.linalg.norm(pB) <= Delta:
            p = pB
        elif np.linalg.norm(pU) >= Delta:
            p = Delta * pU / np.linalg.norm(pU)
        else:
            p_diff = pB - pU
            a = np.dot(p_diff, p_diff)
            b_ = 2 * np.dot(pU, p_diff)
            c = np.dot(pU, pU) - Delta**2
            tau = (-b_ + np.sqrt(b_**2 - 4*a*c)) / (2*a)
            p = pU + tau * (pB - pU)
        
        ared = f(x) - f(x + p)
        pred = -(g.T @ p + 0.5 * p.T @ Q @ p)
        rho = ared / pred
        
        if rho < 0.25:
            Delta *= 0.25
        elif rho > 0.75 and np.linalg.norm(p) >= 0.99 * Delta:
            Delta = min(2 * Delta, Delta_max)
        
        if rho > eta:
            x = x + p
    
    return hist, x

# Plotting Helpers
def plot_function_value(f_sd, f_tr, ylabel, title, fname_prefix):
    plt.figure()
    plt.plot(f_sd, label="Steepest Descent")
    plt.plot(f_tr, label="Trust Region (Dogleg)")
    plt.xlabel("Iteration")
    plt.ylabel(ylabel)
    plt.title(title)
    plt.legend()
    plt.grid(True)
    plt.savefig(f"{fname_prefix}_linear.png", dpi=300)


def plot_grad_norm(g_sd, g_tr, ylabel, title, fname_prefix):
    # Linear scale
    plt.figure()
    plt.plot(g_sd, label="Steepest Descent")
    plt.plot(g_tr, label="Trust Region (Dogleg)")
    plt.xlabel("Iteration")
    plt.ylabel(ylabel)
    plt.title(f"{title} (Linear Scale)")
    plt.legend()
    plt.grid(True)
    plt.savefig(f"{fname_prefix}_linear.png", dpi=300)

    # Log scale
    plt.figure()
    plt.plot(g_sd, label="Steepest Descent")
    plt.plot(g_tr, label="Trust Region (Dogleg)")
    plt.xlabel("Iteration")
    plt.ylabel(ylabel)
    plt.title(f"{title} (Log Scale)")
    plt.yscale('log')
    plt.ylim(1e-20, 1e1)
    plt.legend()
    plt.grid(True, which="both")
    plt.savefig(f"{fname_prefix}_log.png", dpi=300)

# Main Execution
if __name__ == "__main__":
    print("=" * 80)
    print("--- Running Script 1: TR (Dogleg) vs Steepest Descent ---")
    print("=" * 80)

    # Run algorithms
    print("\nRunning Trust Region (Delta_0 = 1.0) to match assignment setup...")
    tr_hist, x_tr = trust_region_dogleg(x0, Delta0=1.0, Delta_max=100.0, tol=1e-8, maxiter=200)
    print(f"Trust Region converged at iteration {max(len(tr_hist['f'])-1,1)}")

    print("\nRunning Steepest Descent (max_iter = 200)...")
    sd_hist, x_sd = steepest_descent(x0, tol=1e-8, maxiter=200)
    if len(sd_hist['f']) == 200:
        print("Steepest Descent reached max_iter (200) without converging.")
    else:
        print(f"Steepest Descent converged at iteration {len(sd_hist['f'])}.")

    # ------------------------------------------------------
    # Summary Output
    # ------------------------------------------------------
    print("-" * 80)
    print(f"Trust Region completed in {max(len(tr_hist['f'])-1,1)} iterations.")
    print(f"Steepest Descent completed in {len(sd_hist['f'])} iterations.\n")

    print(f"Final Trust Region x* =\n {np.round(x_tr, 6)}")
    print(f"Final Steepest Descent x* =\n {np.round(x_sd, 6)}\n")

    print(f"Final function value (TR) = {tr_hist['f'][-1]:.6f}")
    print(f"Final function value (SD) = {sd_hist['f'][-1]:.6f}")
    print(f"||grad(x_TR)|| = {tr_hist['g'][-1]:.3e}")
    print(f"||grad(x_SD)|| = {sd_hist['g'][-1]:.3e}")
    print(f"Difference between x_TR and x_SD = {np.linalg.norm(x_tr - x_sd):.3e}")
    print(f"Difference between function values = {abs(tr_hist['f'][-1] - sd_hist['f'][-1]):.3e}")
    print("-" * 80)

    # Plots
    plot_function_value(sd_hist['f'], tr_hist['f'],
                        "Function Value", "Function Value vs Iteration", "function_value")
    
    plot_grad_norm(sd_hist['g'], tr_hist['g'],
                   "||Gradient||", "Gradient Norm vs Iteration", "grad_norm")

    plt.figure()
    plt.plot(tr_hist['Delta'], label="Trust Region Radius Δ")
    plt.xlabel("Iteration")
    plt.ylabel("Δ (Trust Region Radius)")
    plt.title("Trust Region Radius vs Iteration")
    plt.legend()
    plt.grid(True)
    plt.savefig("delta_vs_iteration.png", dpi=300)