project2
unknown
python
9 months ago
9.3 kB
14
Indexable
import numpy as np
import matplotlib.pyplot as plt
### Function Definitions ###
# Function 1:
def function1(x):
# Sum (i x_i^2) from i=1 to 100
n = 100
i_vals = np.arange(1, n + 1)
return np.sum(i_vals * x**2)
def dfunction1(x):
# Gradient of function1
n = 100
i_vals = np.arange(1, n + 1)
return 2 * i_vals * x
def hfunction1(x):
# Hessian of function1
n = 100
i_vals = np.arange(1, n + 1)
return np.diag(2 * i_vals)
# Function 2:
b = np.loadtxt('files/fun2_b.txt').reshape(-1,1)
c = np.loadtxt('files/fun2_c.txt').reshape(-1,1)
m = b.shape[0]
n = c.shape[0]
A = np.loadtxt('files/fun2_A.txt').reshape(m,n)
def make_f2(A, b, c):
# A = 500 x 100
# b = 500 x 1
# c = 100 x 1
# x = 100
def f(x):
x = x.reshape(-1,1)
return (c.T@x - np.sum(np.log(b - A@x))).item()
def df(x):
x = x.reshape(-1,1)
return (c + A.T@(1/(b-A@x))).flatten()
def hf(x):
x = x.reshape(-1,1)
r = (b-A@x)**2
m = b.shape[0]
return A.T @ ((1/r)*np.eye(m)) @ A
return f, df, hf # , (100,), (100,100)
function2, dfunction2, hfunction2 = make_f2(A, b, c)
# Function 3:
def function3(x):
# f(x) = 100(x2 - x1^2)^2 + (1 - x1)^2
return 100.0 * (x[1] - x[0]**2)**2 + (1 - x[0])**2
def dfunction3(x):
# Gradient of function3
dfdx1 = -400.0 * x[0] * (x[1] - x[0]**2) - 2.0 * (1 - x[0])
dfdx2 = 200.0 * (x[1] - x[0]**2)
return np.array([dfdx1, dfdx2])
def hfunction3(x):
# Hessian of function3
h11 = 1200.0 * x[0]**2 - 400.0 * x[1] + 2.0
h12 = -400.0 * x[0]
h21 = -400.0 * x[0]
h22 = 200.0
return np.array([[h11, h12], [h21, h22]])
### Algorithm Implementation ###
def backtracking(f, gradf, x, p, alpha, rho=0.5, c=1e-4, max_iter=1000):
# Rho: step size reduction factor
# C: sufficient decrease parameter
# Condition: f(x + alpha * p) <= c * p * gradf(x) * alpha + f(x)
for _ in range(max_iter):
if f(x + alpha * p) <= c * np.dot(p, gradf(x)) * alpha + f(x):
return alpha
alpha *= rho
return alpha # Fallback return
def gradient_descent(x, f, df, alpha = 1.0, tol=1e-6, max_iter=1000):
# Alpha: initial step size/learning rate
# Tolerance for stopping criterion
# Max iterations for safety
xs = [x.copy()]
fs = [f(x)]
for it in range(max_iter):
grad = df(x)
if np.linalg.norm(grad) < tol:
break
p = -grad
alpha = backtracking(f, df, x, p, alpha)
x = x + alpha * p
# Store history
xs.append(x.copy())
fs.append(f(x))
return dict(x=x, xs=xs, fs=fs, nit=it+1)
def newtons_method(x, f, df, hf, E=1e-6, max_steps=1000):
xs = [x.copy()]
fs = [f(x)]
for i in range(max_steps):
grad = df(x)
H = hf(x)
p = -np.linalg.solve(H, grad)
alpha = backtracking(f, df, x, p, 1.0)
x = x + alpha * p
# Store history
xs.append(x.copy())
fs.append(f(x))
if np.linalg.norm(df(x)) < E:
break
# If max_steps reached
return dict(x=x, xs=xs, fs=fs, nit=i+1)
def quasi_newtons_method(x, f, df, E=10e-6, max_steps=1000):
F_0 = np.eye(x.shape[0])
alpha = 1.0
xs, fs = [x.copy()], [f(x)]
for i in range(max_steps):
try:
x = x.squeeze(1)
except:
pass
p = -F_0 @ df(x)
alpha = backtracking(f, df, x, p, alpha)
x1 = x + alpha * p
xs.append(x1.copy())
fs.append(f(x1))
# Update F_0 using BFGS formula
s = (x1 - x).reshape(-1, 1)
y = (df(x1) - df(x)).reshape(-1, 1)
denom = y.T @ s
# BFGS update formula
Fy = F_0 @ y
term1 = ((y.T @ (Fy + s)) / (denom ** 2)) * (s @ s.T)
term2 = -((s @ y.T) @ F_0 + F_0 @ (y @ s.T)) / denom
F_1 = F_0 + term1 + term2
F_0 = F_1
x = x1
if np.linalg.norm(df(x)) < E:
break
return dict(x=x, xs=xs, fs=fs, nit=i+1)
### Plot Function ###
def plot_convergence(results, func_name="Function", save_path=None):
"""
Plot convergence comparison showing only function value vs iterations
Parameters:
results: dict with method names as keys and result dicts as values
Each result dict should have 'fs' (function values) and 'nit' (iterations)
func_name: name of the function being optimized
save_path: optional path to save the figure
"""
plt.figure(figsize=(10, 6))
# Plot function value vs iterations
for method_name, result in results.items():
plt.plot(range(len(result['fs'])), result['fs'],
'o-', linewidth=2, markersize=5,
label=f"{method_name} ({result['nit']} iter)")
plt.xlabel('Iteration', fontsize=14)
plt.ylabel('Function Value', fontsize=14)
plt.title(f'{func_name}: Convergence Comparison', fontsize=16, fontweight='bold')
plt.legend(loc='best', fontsize=12)
plt.grid(True, alpha=0.3, linestyle='--')
plt.tight_layout()
if save_path:
plt.savefig(save_path, dpi=300, bbox_inches='tight')
print(f"Plot saved to {save_path}")
return plt
# import matplotlib.pyplot as plt
def plot_convergence_log(results, func_name="Function", save_path=None, log_scale=True):
"""
Plot convergence comparison showing function value vs iterations.
Parameters:
----------
results : dict
Dictionary with method names as keys and result dicts as values.
Each result dict should have:
- 'fs': list or array of function values
- 'nit': number of iterations
func_name : str, optional
Name of the function being optimized (default is "Function").
save_path : str or None, optional
Path to save the figure. If None, the plot is not saved.
log_scale : bool, optional
If True, use symmetric log scale (symlog) for y-axis to support negative values.
Returns:
-------
plt : matplotlib.pyplot
The plot object.
"""
plt.figure(figsize=(10, 6))
# Plot function value vs iterations
for method_name, result in results.items():
plt.plot(range(len(result['fs'])), result['fs'],
'o-', linewidth=2, markersize=5,
label=f"{method_name} ({result['nit']} iter)")
plt.xlabel('Iteration', fontsize=14)
plt.ylabel('Function Value', fontsize=14)
if log_scale:
plt.yscale('symlog', linthresh=1e-8) # symlog supports negative numbers
plt.title(f'{func_name}: Convergence Comparison', fontsize=16, fontweight='bold')
plt.legend(loc='best', fontsize=12)
plt.grid(True, alpha=0.3, linestyle='--')
plt.tight_layout()
if save_path:
plt.savefig(save_path, dpi=300, bbox_inches='tight')
print(f"Plot saved to {save_path}")
return plt
### Generate initial point ###
def generate_initial_point(n):
return np.random.randn(n)
def generate_random_feasible_simple(A, b, c, num_tries=10):
n = c.shape[0]
for attempt in range(num_tries):
# Scale: 1, 0.5, 0.25, ...
scale = 1 / (2 ** attempt)
x_trial = np.random.randn(n) * scale
residual = b - A @ x_trial.reshape(-1, 1)
if np.all(residual > 0):
return x_trial
return np.zeros(n)
def main():
# Set random seed for reproducibility
np.random.seed(42)
x0 = generate_initial_point(100)
results = {
"Gradient Descent": gradient_descent(x0, function1, dfunction1),
"Newton's Method": newtons_method(x0, function1, dfunction1, hfunction1),
"Quasi-Newton's Method": quasi_newtons_method(x0, function1, dfunction1)
}
plot_convergence_log(results, func_name="Function 1", save_path="function1_convergence.png")
plt.show()
print("Predicted f1 results:")
for method, result in results.items():
print(f"{method}: {function1(result['x'])}")
print("\n")
x0 = generate_random_feasible_simple(A, b, c)
results = {
"Gradient Descent": gradient_descent(x0, function2, dfunction2),
"Newton's Method": newtons_method(x0, function2, dfunction2, hfunction2),
"Quasi-Newton's Method": quasi_newtons_method(x0, function2, dfunction2)
}
plot_convergence(results, func_name="Function 2", save_path="function2_convergence.png")
plt.show()
print("Predicted f2 Results:")
for method, result in results.items():
print(f"{method}: {function2(result['x'])}")
print("\n")
x0 = generate_initial_point(2)
results = {
"Gradient Descent": gradient_descent(x0, function3, dfunction3),
"Newton's Method": newtons_method(x0, function3, dfunction3, hfunction3),
"Quasi-Newton's Method": quasi_newtons_method(x0, function3, dfunction3)
}
plot_convergence_log(results, func_name="Function 3", save_path="function3_convergence.png")
plt.show()
print("Predicted f3 Results:")
for method, result in results.items():
print(f"{method}: {function3(result['x'])}")
print("\n")
if __name__ == "__main__":
main()Editor is loading...
Leave a Comment