Code_Problem_Set_4_Math.py

import numpy as np
import matplotlib.pyplot as plt
import time

# Constants
c = 10
n = 2

# x1 and x2 vectors
x1_vals = np.linspace(-1, 1, 100)
x2_vals = np.linspace(-1, 1, 100)

# grid of (x1, x2) points
X1, X2 = np.meshgrid(x1_vals, x2_vals)

# Diagonal entries for C
# For i=1: c^((1-1)/(2-1)) = 10^0 = 1
# For i=2: c^((2-1)/(2-1)) = 10^1 = 10
diag_vals = [c**((i-1)/(n-1)) for i in range(1, n + 1)]

# Create the 2x2 C matrix
C = np.diag(diag_vals)

# fsq(x) = x^T C x
xt_C_x = C[0, 0] * X1**2 + C[1, 1] * X2**2

Z_sq = xt_C_x
# fhole(x) = 1 - exp(x^T C x)
Z_hole = 1 - np.exp(xt_C_x)



fig = plt.figure(figsize=(14, 7))

#Subplots
ax1 = fig.add_subplot(1, 2, 1, projection='3d')
ax1.plot_surface(X1, X2, Z_sq, cmap='viridis')
ax1.set_title('$f_{sq}(x) = x^T C x$')
ax1.set_xlabel('$x_1$')
ax1.set_ylabel('$x_2$')
ax1.set_zlabel('$f(x)$')


ax2 = fig.add_subplot(1, 2, 2, projection='3d')
ax2.plot_surface(X1, X2, Z_hole, cmap='plasma')
ax2.set_title('$f_{hole}(x) = 1 - exp(x^T C x)$')
ax2.set_xlabel('$x_1$')
ax2.set_ylabel('$x_2$')
ax2.set_zlabel('$f(x)$')

plt.show()


#Gradient Formulas
#Gradient is ∇f = 2xᵀC 
def grad_f_sq(x, C):
    #transpose of c is itself since it's diagonal
    return 2 * C.T @ x # @ matrix multiplication

# --- Define the gradient for f_hole ---
# gradient is ∇f = -2xᵀC * exp(xᵀCx).
# As a column vector, this is ∇f = -2Cx * exp(xᵀCx).
def grad_f_hole(x, C):
    xt_C_x = x.T @ C @ x
    exp_term = np.exp(xt_C_x)
    return -2 * C @ x * exp_term

# Gradient Descent
def run_gradient_descent(grad_func, title, C, alpha=0.01, n_iterations=10):
    x = np.array([[1.0], [1.0]])
    
    print(f"--- {title} (Gradient Descent) ---")
    print(f"Step 0: x = {x.flatten()}")
    
    start_time = time.time()
    #print("Start_time:", start_time)
    for i in range(n_iterations):
        grad = grad_func(x, C)
        x = x - alpha * grad
        print(f"Step {i+1}: x = {x.flatten()}")
    end_time = time.time()
    
    print(f"Total time: {end_time - start_time:.6f} seconds\n")

print("grad_f_sq =", run_gradient_descent(grad_f_sq, "f_sq Gradient Descent", C))
print("grad_f_hole =", run_gradient_descent(grad_f_hole, "f_hole Gradient Descent", C))
# Hessians
def hess_f_sq(x, C):
    return 2 * C.T

def hess_f_hole(x, C):
    xt_C_x = float(x.T @ C @ x)
    exp_term = np.exp(xt_C_x)
    xt_C = x.T @ C
    # This is the outer product of the vector Cx with its transpose
    outer_product = xt_C @ xt_C
    return -4 * exp_term * outer_product



def run_newtons_method(grad_func, hess_func, title, C, n_iterations=10):
    x = np.array([[1.0], [1.0]])
    
    print(f"--- {title} (Newton's Method) ---")
    print(f"Step 0: x = {x.flatten()}")
    
    start_time = time.time()
    for i in range(n_iterations):
        grad = grad_func(x, C)
        hess = hess_func(x, C)
        # Newton's update rule: x_k+1 = x_k - H⁻¹ * ∇f
        # We solve the linear system H * Δx = ∇f instead of inverting H
        delta_x = np.linalg.solve(hess, grad)
        x = x - delta_x
        
        print(f"Step {i+1}: x = {x.flatten()}")
    end_time = time.time()
    
    print(f"Total time: {end_time - start_time:.6f} seconds\n")


## Run and Compare
run_gradient_descent(grad_f_sq, "f_sq", C)
run_newtons_method(grad_f_sq, hess_f_sq, "f_sq", C)

run_gradient_descent(grad_f_hole, "f_hole", C, alpha=0.01)
run_newtons_method(grad_f_hole, hess_f_hole, "f_hole", C)