ND/1D Optimization

GoldenSection.py

import numpy as np
import matplotlib.pyplot as plt

# TODO: Change us!
f = lambda x: (x - 2.2) ** 2
a = -3
b = 8


brackets = []
gs = (np.sqrt(5) - 1) / 2
m1 = a + (1 - gs) * (b - a)
m2 = a + gs * (b - a)
f1 = f(m1)
f2 = f(m2)
count = 0
# Begin your modifications below here

while b - a >= 1e-5:
    brackets.append([a, m1, m2, b])
    if f1 > f2:
        a = m1
        m1 = m2
        f1 = f2
        h1 = b - a
        m2 = a + gs * h1
        f2 = f(m2)
    if f1 < f2:
        b = m2
        m2 = m1
        f2 = f1
        m1 = a + (1 - gs) * (b - a)
        f1 = f(m1)
    print(f"Step: {count} | L: {a} | R: {b}")
    count += 1

# End your modifications above here

# Plotting code below, no need to modify
x = np.linspace(-10, 10)
plt.plot(x, f(x))

brackets = np.array(brackets)
names=['a', 'm1', 'm2', 'b']
for i in range(4):
    plt.plot(brackets[:, i], 3 * np.arange(len(brackets)), 'o-', label=names[i])
plt.legend()

NDOptimize.py

import argparse
import numpy as np
from sympy.abc import x, y
from sympy.parsing.sympy_parser import parse_expr
from sympy import N, Symbol, Matrix, Derivative, linsolve, lambdify


parser = argparse.ArgumentParser(description="Find the gradient and Hessian matrix for a function\
                                 from R2 to R")

parser.add_argument('-f', dest='f', type=str, nargs='?',
                    help='a symbolic expression for the function f(x, y)')
parser.add_argument('-x', dest='x', type=float, nargs='?',
                    help='function input x')
parser.add_argument('-y', dest='y', type=float, nargs='?',
                    help='function input y')
args = parser.parse_args()

# Symbolic computations
f = parse_expr(args.f)
grad = Matrix([Derivative(f, x, evaluate=True).simplify(),
               Derivative(f, y, evaluate=True).simplify()])
H = grad.jacobian([x, y])
print(f"Gradient\n{grad}")
print(f"Hessian\n{H}")

# Numerical computations
print("Function f")
print(np.array(f.subs([(x, args.x), (y, args.y)])))
print("Gradient")
print(np.array(grad.subs([(x, args.x), (y, args.y)])))
print("Hessian Matrix")
print(np.array(H.subs([(x, args.x), (y, args.y)])))
print("Newton Step")
A = np.array(H.subs([(x, args.x), (y, args.y)]).tolist(), dtype=np.float64)
b = np.array(grad.subs([(x, args.x), (y, args.y)]).tolist(), dtype=np.float64)
print(np.array([args.x, args.y]).reshape(-1, 1) + np.linalg.solve(A, -b))
print("Steepest Descent")
print(np.array(-(grad.subs([(x, args.x), (y, args.y)]))))

Newton1D.py

import numpy as np

# complete the function below
def dfunc(x):
    # Add your code here
    return np.exp(-x ** 2) * (-1 + 2 * x ** 2 - np.cos(x) + 2 * x * np.sin(x))

# complete the function below
def d2func(x):
    # Add your code here
    return np.exp(-x ** 2) * (6 * x - 4 * x ** 3 + 4 * x * np.cos(x) + np.sin(x) * (3 - 4 * x **2))

# run Newton's Method
xk = x0 - dfunc(x0) / d2func(x0)
guess = [x0, xk]
while np.abs(dfunc(xk)) > tol:
    xk -= dfunc(xk) / d2func(xk)
    guess.append(xk)

newton_guesses = np.array(guess)

SteepestDescent.py

import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt

def f(r):
    x, y = r
    return 3 +((x**2)/8) + ((y**2)/8) - np.sin(x)*np.cos((2**-0.5)*y)

def obj(alpha, r, s):
    return f(r + alpha * s)

def grad(r):
    x, y = r
    gradx = 1 / 4 * (x - 4 * np.cos(x) * np.cos(y / np.sqrt(2)))
    grady = 1 / 4 * (y + 2 * np.sqrt(2) * np.sin(x) * np.sin(y / np.sqrt(2)))
    return np.array([gradx, grady])

def H(r):
    x, y = r
    return np.array([
        [1 / 4 + np.cos(y / np.sqrt(2)) * np.sin(x), np.cos(x) * np.sin(y / np.sqrt(2)) / np.sqrt(2)],
        [np.cos(x) * np.sin(y / np.sqrt(2)) / np.sqrt(2), 1 / 4 + (1 / 2) * np.cos(y / np.sqrt(2)) * np.sin(x)]
    ])

r1 = [[r_init[0], r_init[1]]]
r_sd = r_init.copy()
iteration_count_sd = 0
while np.linalg.norm(grad(r_sd), 2) >= stop:
    iteration_count_sd += 1
    s = -grad(r_sd)
    alpha = opt.minimize_scalar(obj, args=(r_sd, s)).x
    r_sd += alpha * s
    r1.append([r_sd[0], r_sd[1]])

r2 = [[r_init[0], r_init[1]]]
r_newton = r_init.copy()
iteration_count_newton = 0
while np.linalg.norm(grad(r_newton), 2) >= stop:
    iteration_count_newton += 1
    s = np.linalg.solve(H(r_newton), -grad(r_newton))
    r_newton += s
    r2.append([r_newton[0], r_newton[1]])

r1 = np.array(r1) - r_sd
r2 = np.array(r2) - r_newton

plt.plot(np.arange(len(r1)), np.log(np.linalg.norm(r1, axis=1)))
plt.plot(np.arange(len(r2)), np.log(np.linalg.norm(r2, axis=1)))
plt.xlabel("Number of Iterations")
plt.ylabel("Error")
plt.title("Error vs # of Iterations")
plt.legend()