diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
new file mode 100644
index 0000000..96e3141
--- /dev/null
+++ b/.github/workflows/test.yml
@@ -0,0 +1,17 @@
+name: test
+
+on:
+ push:
+ pull_request:
+
+jobs:
+ pytest:
+ runs-on: ubuntu-latest
+ steps:
+ - uses: actions/checkout@v4
+ - uses: actions/setup-python@v5
+ with:
+ python-version: "3.9"
+ cache: pip
+ - run: python -m pip install -r requirements-dev.txt
+ - run: MPLBACKEND=Agg pytest -q
diff --git a/.gitignore b/.gitignore
index d2d6f36..d5abaa0 100644
--- a/.gitignore
+++ b/.gitignore
@@ -33,3 +33,8 @@ nosetests.xml
.mr.developer.cfg
.project
.pydevproject
+
+# Local development
+.venv/
+__pycache__/
+.pytest_cache/
diff --git a/README.md b/README.md
index 8f9ffd2..81e2d3f 100644
--- a/README.md
+++ b/README.md
@@ -1,8 +1,23 @@
-gradient
-========
-
-optimisation practices
-- Newton, modified Newton, Quasi Newton (rank one with backtracking)
-
-
-
+# gradient
+
+Optimization path comparisons:
+
+- Newton and modified Newton
+- symmetric-rank-one quasi-Newton
+- BFGS
+- Fletcher–Reeves nonlinear conjugate gradient
+- backtracking and interpolating Armijo line searches
+
+The original figure layout and plotting style are preserved. Running `main.py`
+regenerates the three figures with the corrected numerical paths.
+
+```sh
+python3 -m venv .venv
+.venv/bin/pip install -r requirements-dev.txt
+.venv/bin/pytest
+MPLBACKEND=Agg .venv/bin/python main.py
+```
+
+
+
+
diff --git a/cg.png b/cg.png
index 20004f4..66ffb0e 100644
Binary files a/cg.png and b/cg.png differ
diff --git a/drawFunc.py b/drawFunc.py
index 90ff8a6..1a96ecb 100644
--- a/drawFunc.py
+++ b/drawFunc.py
@@ -1,21 +1,24 @@
-import numpy as np
+"""Original contour plotting helpers."""
+
import matplotlib.pyplot as plt
-from func import *
+import numpy as np
-def draw(f_x):
-# allocate grids
- dx = np.arange(-3.0, 2.5, .02)
- dy = np.arange(-4.5, 4.5, .02)
- X, Y = np.meshgrid(dx, dy)
-# eval function f_x
- Z = [f_x(np.matrix([[x1], [x2]])) \
- for (x1, x2) in zip(np.hstack(X), np.hstack(Y))]
- Z = np.array(Z)
- Z = Z.reshape(np.shape(X))
+def draw(
+ f_x,
+ filled=True,
+ x_bounds=(-3.0, 2.5),
+ y_bounds=(-4.5, 4.5),
+ levels=None,
+):
+ x_values = np.arange(x_bounds[0], x_bounds[1], 0.02)
+ y_values = np.arange(y_bounds[0], y_bounds[1], 0.02)
+ X, Y = np.meshgrid(x_values, y_values)
+ points = np.column_stack((X.ravel(), Y.ravel()))
+ Z = np.fromiter((f_x(point) for point in points), dtype=float).reshape(X.shape)
-# take this pencil
- levels = [0.0, 1, 2, 5, 20, 30, 50, 80, 100, 200, 400]
- #cs = plt.contour(X, Y, Z, levels, colors='k')
- cs = plt.contourf(X, Y, Z, 20, cmap=plt.cm.YlGnBu)
- #plt.clabel(cs, fmt='%.1f', inline=1)
+ if filled:
+ return plt.contourf(X, Y, Z, 20, cmap=plt.cm.YlGnBu)
+ contours = plt.contour(X, Y, Z, levels=levels, colors="black")
+ plt.clabel(contours, fmt="%.3f", inline=True)
+ return contours
diff --git a/func.py b/func.py
index 6a32d9e..01293ae 100644
--- a/func.py
+++ b/func.py
@@ -1,59 +1,64 @@
+"""Objective functions used by the optimisation examples."""
+
import numpy as np
-class quadratic:
-# trivial (solving linear system)
- A = np.zeros([2, 2])
- b = np.zeros([2, 1])
+def _vector(x):
+ """Return *x* as a one-dimensional floating-point array."""
+ return np.asarray(x, dtype=float).reshape(-1)
- def __init__(self):
- pass
+class quadratic:
def __init__(self, A, b):
- self.A = A
- self.b = b
+ self.A = np.asarray(A, dtype=float)
+ self.b = _vector(b)
+ if self.A.shape != (self.b.size, self.b.size):
+ raise ValueError("A must be square and match the size of b")
+ if not np.all(np.isfinite(self.A)) or not np.all(np.isfinite(self.b)):
+ raise ValueError("A and b must contain only finite values")
+ if not np.allclose(self.A, self.A.T):
+ raise ValueError("A must be symmetric so A @ x - b is the gradient")
def f_x(self, x):
- f = 0.5 * np.dot(np.dot(x.T, self.A), x) - np.dot(self.b.T, x)
- return f[0, 0]
+ x = _vector(x)
+ return float(0.5 * x @ self.A @ x - self.b @ x)
def g_x(self, x):
- return np.dot(self.A, x) - self.b
+ return self.A @ _vector(x) - self.b
def G_x(self, x):
- return np.asmatrix(self.A)
+ return self.A.copy()
+
+ def about_alpha(self, x, s):
+ return lambda alpha: self.f_x(_vector(x) + alpha * _vector(s))
+
+ def about_alpha_prime(self, x, s):
+ return lambda alpha: float(self.g_x(_vector(x) + alpha * _vector(s)) @ _vector(s))
class rosenbrock:
-
-# f(x_0, x_1) = 10 * (x_1 - x_0^2)^2 + (x_0 - 1)^2
- def __init__(self):
- pass
+ """The two-dimensional Rosenbrock function with coefficient 10."""
def f_x(self, x):
- f = 10.0 * (x[1] - (x[0])**2)**2 + (1-x[0])**2
- return f[0, 0]
+ x0, x1 = _vector(x)
+ return float(10.0 * (x1 - x0**2) ** 2 + (1.0 - x0) ** 2)
def g_x(self, x):
- g0 = -40.0 * (x[0]*x[1] - x[0]**3) - 2.0 + 2.0*x[0]
- g1 = 20.0 * (x[1] - x[0]**2)
- return np.matrix([[g0[0,0]], [g1[0,0]]])
+ x0, x1 = _vector(x)
+ return np.array([
+ -40.0 * x0 * (x1 - x0**2) + 2.0 * (x0 - 1.0),
+ 20.0 * (x1 - x0**2),
+ ])
def G_x(self, x):
- G00 = 120.0 * x[0]**2 - 40.0 * x[1] + 2.0
- G01 = -40.0*x[0]
- G10 = -40.0*x[0]
- G11 = 20.0
- return np.matrix([[G00[0,0], G01[0,0]], [G10[0,0], G11]])
+ x0, x1 = _vector(x)
+ return np.array([
+ [120.0 * x0**2 - 40.0 * x1 + 2.0, -40.0 * x0],
+ [-40.0 * x0, 20.0],
+ ])
def about_alpha(self, x, s):
- # f(a) = f(x + a*s) at point x, direction s
- def along_s(alpha):
- return self.f_x(x+s*alpha)
- return along_s
+ return lambda alpha: self.f_x(_vector(x) + alpha * _vector(s))
def about_alpha_prime(self, x, s):
- def along_s(alpha):
- p = np.dot(self.g_x(x+alpha*s).T, s)
- return p[0, 0]
- return along_s
+ return lambda alpha: float(self.g_x(_vector(x) + alpha * _vector(s)) @ _vector(s))
diff --git a/main.py b/main.py
index 1685be6..4a060b1 100644
--- a/main.py
+++ b/main.py
@@ -1,100 +1,147 @@
-import optMethods as opt
+"""Regenerate the repository's three original-style path figures."""
+
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+
import drawFunc
import func
-import numpy as np
-import matplotlib.pyplot as plt
-import pdb
-
-def quadratic():
-
- A = np.matrix([
- [3., 2.],
- [2., 6.]])
- b = np.matrix([[1.], [-5.]])
- obj = func.quadratic(A,b)
-
- plt.figure()
- drawFunc.draw(obj.f_x)
-
- x = np.matrix([[0], [2]])
- track = opt.newton(x, obj, 0, 1) # just one iteration
- p1, = plt.plot(track[0,:].tolist()[0], track[1,:].tolist()[0])
-
- track = opt.quasi_newton(x, obj, 100)
- p2, = plt.plot(track[0,:].tolist()[0], track[1,:].tolist()[0],
- 'o-', linewidth=2.0)
-
- plt.legend([p1, p2], ['Newton', 'Quasi_Newton_Rank_One'])
-
- plt.show()
-
-def rosen():
- obj = func.rosenbrock()
-
- plt.figure()
-
- drawFunc.draw(obj.f_x)
- start_point = np.matrix([[-1.5], [-4.0]])
- plt.annotate('Start', xy=(-1.5, -4.0), xytext=(-1.8, -3.5),
- arrowprops=dict(facecolor='black', shrink=0.02, frac=0.5))
- plt.annotate('Optimal', xy=(1, 1), xytext=(1.1, 1.4),
- arrowprops=dict(facecolor='black', shrink=0.02, frac=0.5))
-
- #v = 0.0
- #str1 = "Newton"
- #track = opt.newton(start_point, obj, v)
- #p1, = plt.plot(track[0,:].tolist()[0], track[1,:].tolist()[0],
- # 'o-', linewidth=2.0)
-
- v = 0.1
- str2 = "Modified_Newton %.2f"%v
- track = opt.newton(start_point, obj, v)
- p2, = plt.plot(track[0,:].tolist()[0], track[1,:].tolist()[0],
- 'o-', linewidth=2.0)
-
- v = 1.5
- str3 = "Modified_Newton %.2f"%v
- track = opt.newton(start_point, obj, v, 50)
- p3, = plt.plot(track[0,:].tolist()[0], track[1,:].tolist()[0],
- 'o-', linewidth=2.0)
-
- str4 = "Quasi_Newton_trivial"
- track = opt.quasi_newton(start_point, obj, 150, 1)
- p4, = plt.plot(track[0,:].tolist()[0], track[1,:].tolist()[0],
- 'o-', linewidth=2.0)
-
- str5 = "Quasi_Newton_Armijo"
- track = opt.quasi_newton(start_point, obj, 450, 0)
- p5, = plt.plot(track[0,:].tolist()[0], track[1,:].tolist()[0],
- 'o-', linewidth=2.0)
-
- str6 = "BFGS_trivial"
- track = opt.BFGS(start_point, obj, 350, 1)
- p6, = plt.plot(track[0,:].tolist()[0], track[1,:].tolist()[0],
- 'o-', linewidth=2.0)
-
- str7 = "BFGS_Armijo"
- track = opt.BFGS(start_point, obj, 350, 0)
- p7, = plt.plot(track[0,:].tolist()[0], track[1,:].tolist()[0],
- 'o-', linewidth=2.0)
-
- v = 0.001
- str8 = "CG_FR_%.5f"%v
- track = opt.fletcher_reeves(start_point, obj, iteration=30, alpha=v)
- p8, = plt.plot(track[0,:].tolist()[0], track[1,:].tolist()[0],
- 'o-', linewidth=2.0)
-
- v = 0.01
- str9 = "CG_FR_%.5f"%v
- track = opt.fletcher_reeves(start_point, obj, iteration=30, alpha=v)
- p9, = plt.plot(track[0,:].tolist()[0], track[1,:].tolist()[0],
- 'o-', linewidth=2.0)
- #plt.legend([p1, p2, p3, p4, p5, p6, p7],\
- # [str1, str2, str3, str4, str5, str6, str7])
-
- plt.legend([p2, p3, p4, p5, p6, p7, p8, p9],\
- [str2, str3, str4, str5, str6, str7, str8, str9])
- plt.show()
-
-#quadratic()
-rosen()
+import optMethods as opt
+
+
+ROOT = Path(__file__).resolve().parent
+
+
+def _plot(track, *args, **kwargs):
+ return plt.plot(track[0], track[1], *args, **kwargs)[0]
+
+
+def quadratic(output=ROOT / "quadratic.png"):
+ objective = func.quadratic([[3.0, 2.0], [2.0, 6.0]], [1.0, -5.0])
+ start = np.array([0.0, 2.0])
+
+ plt.figure(figsize=(8, 6))
+ drawFunc.draw(
+ objective.f_x,
+ filled=False,
+ x_bounds=(-1.5, 2.3),
+ y_bounds=(-3.0, 3.5),
+ levels=np.arange(0, 68, 4),
+ )
+ newton_path = opt.newton(start, objective, iterations=1)
+ sr1_path = opt.quasi_newton(start, objective, iteration=100)
+ newton_line = _plot(newton_path, "b-")
+ sr1_line = _plot(
+ sr1_path,
+ "go-",
+ linewidth=2.0,
+ markeredgecolor="black",
+ markeredgewidth=0.5,
+ )
+ plt.legend(
+ [newton_line, sr1_line],
+ ["Newton", "Quasi_Newton_Rank_One"],
+ loc="upper right",
+ )
+ plt.xlim(-1.5, 2.3)
+ plt.ylim(-3.0, 3.5)
+ plt.tight_layout()
+ plt.savefig(output, dpi=100)
+ plt.close()
+
+
+def quasi_newton(output=ROOT / "quasi_newton.png"):
+ objective = func.rosenbrock()
+ start = np.array([-1.0, 2.0])
+
+ plt.figure(figsize=(8, 6))
+ drawFunc.draw(
+ objective.f_x,
+ filled=False,
+ x_bounds=(-1.5, 2.0),
+ y_bounds=(-3.0, 3.5),
+ levels=[0, 10, 25, 50, 100, 150, 200, 250],
+ )
+ plt.annotate(
+ "Start",
+ xy=start,
+ xytext=(-1.2, 2.2),
+ arrowprops=dict(facecolor="black", shrink=0.02),
+ )
+ plt.annotate(
+ "Optimal",
+ xy=(1, 1),
+ xytext=(0.5, 1.2),
+ arrowprops=dict(facecolor="black", shrink=0.02),
+ )
+
+ paths = [
+ (opt.newton(start, objective, iterations=20), "bo-", "Newton"),
+ (opt.newton(start, objective, modified=0.1, iterations=80), "go-", "Modified_Newton 0.10"),
+ (opt.newton(start, objective, modified=1.5, iterations=100), "ro-", "Modified_Newton 1.50"),
+ (opt.quasi_newton(start, objective, iteration=250), "co-", "Quasi_Newton_Rank_One"),
+ ]
+ lines = [
+ _plot(path, style, linewidth=2.0, markersize=4, markeredgecolor="black", markeredgewidth=0.4)
+ for path, style, _ in paths
+ ]
+ plt.legend(lines, [label for _, _, label in paths], loc="upper right")
+ plt.xlim(-1.5, 2.0)
+ plt.ylim(-3.0, 3.5)
+ plt.tight_layout()
+ plt.savefig(output, dpi=100)
+ plt.close()
+
+
+def comparison(output=ROOT / "cg.png"):
+ objective = func.rosenbrock()
+ start = np.array([-1.5, -4.0])
+
+ plt.figure(figsize=(10, 8.5))
+ drawFunc.draw(objective.f_x, filled=True)
+ plt.annotate(
+ "Start",
+ xy=start,
+ xytext=(-1.8, -3.5),
+ arrowprops=dict(facecolor="black", shrink=0.02),
+ )
+ plt.annotate(
+ "Optimal",
+ xy=(1, 1),
+ xytext=(1.1, 1.4),
+ arrowprops=dict(facecolor="black", shrink=0.02),
+ )
+
+ paths = [
+ (opt.newton(start, objective, modified=0.1, iterations=80), "Modified_Newton 0.10"),
+ (opt.newton(start, objective, modified=1.5, iterations=100), "Modified_Newton 1.50"),
+ (opt.quasi_newton(start, objective, iteration=250, interpolation=1), "Quasi_Newton_trivial"),
+ (opt.quasi_newton(start, objective, iteration=450), "Quasi_Newton_Armijo"),
+ (opt.BFGS(start, objective, iteration=350, interpolation=1), "BFGS_trivial"),
+ (opt.BFGS(start, objective, iteration=350), "BFGS_Armijo"),
+ (opt.fletcher_reeves(start, objective, iteration=80, alpha=0.001), "CG_FR_0.00100"),
+ (opt.fletcher_reeves(start, objective, iteration=80, alpha=0.01), "CG_FR_0.01000"),
+ ]
+ lines = [
+ _plot(path, "o-", linewidth=2.0, markersize=4, markeredgecolor="black", markeredgewidth=0.4)
+ for path, _ in paths
+ ]
+ plt.legend(lines, [label for _, label in paths], loc="upper right")
+ plt.xlim(-3.0, 2.5)
+ plt.ylim(-4.5, 4.5)
+ plt.tight_layout()
+ plt.savefig(output, dpi=120)
+ plt.close()
+
+
+def generate_figures(output_dir=ROOT):
+ output_dir = Path(output_dir)
+ output_dir.mkdir(parents=True, exist_ok=True)
+ quadratic(output_dir / "quadratic.png")
+ quasi_newton(output_dir / "quasi_newton.png")
+ comparison(output_dir / "cg.png")
+
+
+if __name__ == "__main__":
+ generate_figures()
diff --git a/optMethods.py b/optMethods.py
index e5c8380..0ce5bdd 100644
--- a/optMethods.py
+++ b/optMethods.py
@@ -1,179 +1,166 @@
-from numpy import *
-
-# TODO: convergence check
-
-def newton(start_point, obj_fun, modified=0, iterations=5):
-# with second order information
- x = start_point
- track = x
-
- k = 0
- while k < iterations:
-
- if modified > 0:
- vI = modified * matrix([[1, 0], [0, 1]])
- G_ = linalg.inv(obj_fun.G_x(x) + vI) # inverse?
- else:
- G_ = linalg.inv(obj_fun.G_x(x)) # inverse?
- g_ = obj_fun.g_x(x)
- delta = -1.0 * dot(G_, g_)
- x = x + delta
-
- track = concatenate((track, x), axis=1)
- k += 1
-
- return track
-
-def BFGS(start_point, obj_fun, iteration=10, interpolation=0):
-# with first order information
- x = start_point
- track = x
-
- k = 0
- #H = dot(obj_fun.g_x(x), obj_fun.g_x(x).T)
- H = 1.0*eye(2)
- gamma = 1.0
-
- while k < iteration and linalg.norm(gamma) > 1e-10:
-
- p = -dot(H, obj_fun.g_x(x))
-
- if interpolation > 0:
- alpha_k = _backtracking_line_search(obj_fun, x, p)
- else:
- alpha_k = _armijo_line_search(obj_fun, x, p)
-
- s = alpha_k * p
- g_k = obj_fun.g_x(x)
- x = x + alpha_k * p
- g_k_1 = obj_fun.g_x(x)
-
- y = g_k_1 - g_k
-
- z = dot(H, y)
- sTy = dot(s.T, y)
- if sTy > 0:
- H += outer(s, s) * (sTy + dot(y.T, z))[0,0]/(sTy**2) \
- - (outer(z, s) + outer(s, z))/sTy
-
- track = concatenate((track, x), axis=1)
- k += 1
-
- return track
-
-def quasi_newton(start_point, obj_fun, iteration=10, interpolation=0):
-# with first order information
- x = start_point
- track = x
-
- k = 0
- #H = matrix([[.01, 0], [0, .01]])
- H = 1.0*eye(2)
-
- while k < iteration:
-
- s = -1.0 * dot(H, obj_fun.g_x(x))
-
- if interpolation > 0:
- alpha_k = _backtracking_line_search(obj_fun, x, s)
- else:
- alpha_k = _armijo_line_search(obj_fun, x, s)
-
- delta = alpha_k * s
- x_k_1 = x + delta
-
- gamma = obj_fun.g_x(x_k_1) - obj_fun.g_x(x)
- u = delta - dot(H, gamma)
-
- scale_a = dot(u.T, gamma)
- if scale_a == 0: # :(
- scale_a = 0.000001
- H = H + outer(u, u) / scale_a
-
- track = concatenate((track, x), axis=1)
- x = x_k_1
- k += 1
-
- return track
-
-def fletcher_reeves(start_point, obj_fun, iteration=10, alpha=0.1):
- x = start_point
- track = x
-
- k = 0
- while k < iteration:
-
- if k < 1:
- g = obj_fun.g_x(x)
- beta = 0.0
- s = -g
- else:
- g = obj_fun.g_x(x)
-
- beta = dot(g.T, g)/dot(g_old.T, g_old)
- beta = beta[0,0]
-
- s = -g + beta * s
- alpha_k = _armijo_line_search(obj_fun, x, s, alpha0=alpha)
- g_old = obj_fun.g_x(x)
- x = x + alpha_k * s
-
- track = concatenate((track, x), axis=1)
- k += 1
-
- return track
-
-def _backtracking_line_search(obj_fun, x, s, c=1e-4):
-# dummy
- alpha = 1.0
- p = 0.5 # magic number
- c = 0.0001
- diff = 1
-
- f_alpha = obj_fun.about_alpha(x, s)
- g_alpha = obj_fun.about_alpha_prime(x, s)
- i = 0
- while (diff > 0 and i < 50):
- f_x = obj_fun.f_x(x)
- diff = f_alpha(alpha) - f_alpha(0) - c * alpha * g_alpha(0)
- alpha *= p
- i += 1
- return alpha
-
-
-def _armijo_line_search(obj_fun, x, s, c=1e-4, alpha0=1.0):
-# adapted from scipy.optimisation
- amin = 0.0
- f_alpha = obj_fun.about_alpha(x, s)
- g_alpha = obj_fun.about_alpha_prime(x, s)
-
- if(f_alpha(alpha0) <= f_alpha(0) + c * alpha0 * g_alpha(0)):
- return alpha0
-
- alpha1 = -(g_alpha(0)) * alpha0**2 / \
- 2.0 / (f_alpha(alpha0) - f_alpha(0) - g_alpha(0) * alpha0)
- if(f_alpha(alpha1) <= f_alpha(0) + c * alpha1 * g_alpha(0)):
- return alpha1
-
- while alpha1 > amin:
-
- factor = alpha0**2 * alpha1**2 * (alpha1 - alpha0)
- a = alpha0**2 * (f_alpha(alpha1) - f_alpha(0) - g_alpha(0) * alpha1) - \
- alpha1**2 * (f_alpha(alpha0) - f_alpha(0) - g_alpha(0) * alpha0)
- a = a / factor
-
- b = -alpha0**3 * (f_alpha(alpha1) - f_alpha(0) - g_alpha(0) * alpha1)+ \
- alpha1**3 * (f_alpha(alpha0) - f_alpha(0) - g_alpha(0) * alpha0)
- b = b / factor
-
- alpha2 = (-b + sqrt(abs(b**2 - 3 * a * g_alpha(0)))) / (3.0 * a)
- if(f_alpha(alpha2) <= f_alpha(0) + c * alpha2 * g_alpha(0)):
- return alpha2
-
- if(alpha1 - alpha2) > alpha1 / 2.0 or (1 - alpha2/alpha1) < 0.96:
- alpha2 = alpha1 / 2.0
-
- alpha0 = alpha1
- alpha1 = alpha2
-
- return 0.0001
+"""Optimization methods used by the path-comparison figures."""
+
+import numpy as np
+
+
+def _start(start_point):
+ x = np.asarray(start_point, dtype=float).reshape(-1).copy()
+ return x, [x.copy()]
+
+
+def _track(points):
+ return np.column_stack(points)
+
+
+def _gradient(obj_fun, x):
+ return np.asarray(obj_fun.g_x(x), dtype=float).reshape(-1)
+
+
+def newton(start_point, obj_fun, modified=0, iterations=5, tolerance=1e-10):
+ """Newton's method with the tutorial's optional fixed Hessian shift."""
+ x, points = _start(start_point)
+ identity = np.eye(x.size)
+ for _ in range(iterations):
+ gradient = _gradient(obj_fun, x)
+ if np.linalg.norm(gradient) <= tolerance:
+ break
+ hessian = np.asarray(obj_fun.G_x(x), dtype=float) + modified * identity
+ x = x + np.linalg.solve(hessian, -gradient)
+ points.append(x.copy())
+ return _track(points)
+
+
+def BFGS(start_point, obj_fun, iteration=10, interpolation=0, tolerance=1e-10):
+ """Inverse-BFGS with a guarded curvature update."""
+ x, points = _start(start_point)
+ inverse_hessian = np.eye(x.size)
+ gradient = _gradient(obj_fun, x)
+
+ for _ in range(iteration):
+ if np.linalg.norm(gradient) <= tolerance:
+ break
+ direction = -inverse_hessian @ gradient
+ if direction @ gradient >= 0:
+ inverse_hessian = np.eye(x.size)
+ direction = -gradient
+ alpha = _line_search(obj_fun, x, direction, interpolation)
+ step = alpha * direction
+ next_x = x + step
+ next_gradient = _gradient(obj_fun, next_x)
+ y = next_gradient - gradient
+ curvature = float(step @ y)
+ if curvature > 1e-12 * np.linalg.norm(step) * np.linalg.norm(y):
+ rho = 1.0 / curvature
+ identity = np.eye(x.size)
+ left = identity - rho * np.outer(step, y)
+ inverse_hessian = (
+ left @ inverse_hessian @ left.T + rho * np.outer(step, step)
+ )
+ x, gradient = next_x, next_gradient
+ points.append(x.copy())
+ return _track(points)
+
+
+def quasi_newton(
+ start_point, obj_fun, iteration=10, interpolation=0, tolerance=1e-10
+):
+ """Symmetric-rank-one inverse-Hessian quasi-Newton method."""
+ x, points = _start(start_point)
+ inverse_hessian = np.eye(x.size)
+ gradient = _gradient(obj_fun, x)
+
+ for _ in range(iteration):
+ if np.linalg.norm(gradient) <= tolerance:
+ break
+ direction = -inverse_hessian @ gradient
+ if direction @ gradient >= 0:
+ inverse_hessian = np.eye(x.size)
+ direction = -gradient
+ alpha = _line_search(obj_fun, x, direction, interpolation)
+ step = alpha * direction
+ next_x = x + step
+ next_gradient = _gradient(obj_fun, next_x)
+ y = next_gradient - gradient
+ residual = step - inverse_hessian @ y
+ denominator = float(residual @ y)
+ threshold = 1e-8 * np.linalg.norm(residual) * np.linalg.norm(y)
+ if abs(denominator) > threshold:
+ inverse_hessian += np.outer(residual, residual) / denominator
+ x, gradient = next_x, next_gradient
+ # Record the accepted point, not the previous iterate.
+ points.append(x.copy())
+ return _track(points)
+
+
+def fletcher_reeves(
+ start_point, obj_fun, iteration=10, alpha=0.1, tolerance=1e-10
+):
+ """Fletcher-Reeves nonlinear conjugate gradient with safe restarts."""
+ x, points = _start(start_point)
+ gradient = _gradient(obj_fun, x)
+ direction = -gradient
+
+ for _ in range(iteration):
+ if np.linalg.norm(gradient) <= tolerance:
+ break
+ if direction @ gradient >= 0:
+ direction = -gradient
+ step_size = _armijo_line_search(obj_fun, x, direction, alpha0=alpha)
+ next_x = x + step_size * direction
+ next_gradient = _gradient(obj_fun, next_x)
+ denominator = float(gradient @ gradient)
+ beta = float(next_gradient @ next_gradient) / denominator if denominator else 0.0
+ direction = -next_gradient + beta * direction
+ x, gradient = next_x, next_gradient
+ points.append(x.copy())
+ return _track(points)
+
+
+def _line_search(obj_fun, x, direction, interpolation):
+ if interpolation > 0:
+ return _backtracking_line_search(obj_fun, x, direction)
+ return _armijo_line_search(obj_fun, x, direction)
+
+
+def _backtracking_line_search(
+ obj_fun, x, direction, c=1e-4, alpha0=1.0, contraction=0.5
+):
+ """Geometric Armijo backtracking."""
+ if not 0 < c < 1 or not 0 < contraction < 1:
+ raise ValueError("line-search constants must lie between zero and one")
+ alpha = float(alpha0)
+ if not np.isfinite(alpha) or alpha <= 0:
+ raise ValueError("alpha0 must be positive and finite")
+ value = obj_fun.f_x(x)
+ slope = float(_gradient(obj_fun, x) @ direction)
+ if slope >= 0:
+ raise ValueError("line-search direction must be a descent direction")
+ for _ in range(60):
+ if obj_fun.f_x(x + alpha * direction) <= value + c * alpha * slope:
+ # The old implementation contracted once more before returning.
+ return alpha
+ alpha *= contraction
+ raise RuntimeError("backtracking line search failed")
+
+
+def _armijo_line_search(obj_fun, x, direction, c=1e-4, alpha0=1.0):
+ """Safeguarded quadratic-interpolation Armijo search."""
+ if not 0 < c < 1:
+ raise ValueError("c must lie between zero and one")
+ alpha = float(alpha0)
+ if not np.isfinite(alpha) or alpha <= 0:
+ raise ValueError("alpha0 must be positive and finite")
+ value = obj_fun.f_x(x)
+ slope = float(_gradient(obj_fun, x) @ direction)
+ if slope >= 0:
+ raise ValueError("line-search direction must be a descent direction")
+
+ for _ in range(60):
+ trial_value = obj_fun.f_x(x + alpha * direction)
+ if trial_value <= value + c * alpha * slope:
+ return alpha
+ denominator = 2.0 * (trial_value - value - alpha * slope)
+ candidate = -slope * alpha**2 / denominator if denominator > 0 else alpha / 2
+ # Safeguard interpolation so every rejected trial makes real progress.
+ alpha = float(np.clip(candidate, 0.1 * alpha, 0.5 * alpha))
+ raise RuntimeError("Armijo interpolation failed")
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..ae2b0b6
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,4 @@
+[tool.pytest.ini_options]
+pythonpath = ["."]
+testpaths = ["tests"]
+filterwarnings = ["ignore::pyparsing.PyparsingDeprecationWarning"]
diff --git a/quadratic.png b/quadratic.png
index f60320a..59e1822 100644
Binary files a/quadratic.png and b/quadratic.png differ
diff --git a/quasi_newton.png b/quasi_newton.png
index e8fb3d2..594fb33 100644
Binary files a/quasi_newton.png and b/quasi_newton.png differ
diff --git a/requirements-dev.txt b/requirements-dev.txt
new file mode 100644
index 0000000..44dca33
--- /dev/null
+++ b/requirements-dev.txt
@@ -0,0 +1,2 @@
+-r requirements.txt
+pytest==8.4.2
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 0000000..90e52c6
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,2 @@
+matplotlib==3.9.4
+numpy==2.0.2
diff --git a/tests/test_optimizers.py b/tests/test_optimizers.py
new file mode 100644
index 0000000..24a4b59
--- /dev/null
+++ b/tests/test_optimizers.py
@@ -0,0 +1,116 @@
+from pathlib import Path
+
+import matplotlib
+import numpy as np
+import pytest
+from PIL import Image
+
+matplotlib.use("Agg")
+
+import func
+import main
+import optMethods as opt
+
+
+@pytest.fixture
+def quadratic_objective():
+ return func.quadratic([[3.0, 2.0], [2.0, 6.0]], [1.0, -5.0])
+
+
+def test_quadratic_rejects_nonsymmetric_matrix():
+ with pytest.raises(ValueError, match="symmetric"):
+ func.quadratic([[1.0, 2.0], [0.0, 1.0]], [0.0, 0.0])
+
+
+def test_rosenbrock_derivatives_match_finite_differences():
+ objective = func.rosenbrock()
+ x = np.array([-0.7, 1.3])
+ step = 1e-6
+ gradient = np.empty(2)
+ hessian = np.empty((2, 2))
+ for index in range(2):
+ offset = np.zeros(2)
+ offset[index] = step
+ gradient[index] = (objective.f_x(x + offset) - objective.f_x(x - offset)) / (2 * step)
+ hessian[:, index] = (objective.g_x(x + offset) - objective.g_x(x - offset)) / (2 * step)
+ np.testing.assert_allclose(objective.g_x(x), gradient, rtol=1e-5)
+ np.testing.assert_allclose(objective.G_x(x), hessian, rtol=1e-5)
+
+
+def test_newton_reaches_quadratic_optimum_in_one_step(quadratic_objective):
+ path = opt.newton([0, 2], quadratic_objective, iterations=1)
+ expected = np.linalg.solve(quadratic_objective.A, quadratic_objective.b)
+ np.testing.assert_allclose(path[:, -1], expected, atol=1e-12)
+
+
+@pytest.mark.parametrize("method", [opt.BFGS, opt.quasi_newton])
+def test_quasi_newton_methods_converge_on_quadratic(method, quadratic_objective):
+ path = method([0, 2], quadratic_objective, iteration=30)
+ expected = np.linalg.solve(quadratic_objective.A, quadratic_objective.b)
+ np.testing.assert_allclose(path[:, -1], expected, atol=1e-6)
+
+
+def test_sr1_records_the_accepted_iterate(quadratic_objective):
+ path = opt.quasi_newton([0, 2], quadratic_objective, iteration=1)
+ assert path.shape == (2, 2)
+ assert not np.array_equal(path[:, 0], path[:, 1])
+
+
+def test_backtracking_returns_the_step_it_accepted(quadratic_objective):
+ x = np.array([0.0, 2.0])
+ direction = -quadratic_objective.g_x(x)
+ alpha = opt._backtracking_line_search(quadratic_objective, x, direction)
+ slope = quadratic_objective.g_x(x) @ direction
+ assert quadratic_objective.f_x(x + alpha * direction) <= (
+ quadratic_objective.f_x(x) + 1e-4 * alpha * slope
+ )
+ assert quadratic_objective.f_x(x + 2 * alpha * direction) > (
+ quadratic_objective.f_x(x) + 1e-4 * 2 * alpha * slope
+ )
+
+
+def test_interpolating_armijo_returns_an_accepted_step(quadratic_objective):
+ x = np.array([0.0, 2.0])
+ direction = -quadratic_objective.g_x(x)
+ alpha = opt._armijo_line_search(quadratic_objective, x, direction)
+ slope = quadratic_objective.g_x(x) @ direction
+ assert quadratic_objective.f_x(x + alpha * direction) <= (
+ quadratic_objective.f_x(x) + 1e-4 * alpha * slope
+ )
+
+
+def test_line_search_rejects_ascent_direction(quadratic_objective):
+ x = np.array([0.0, 2.0])
+ with pytest.raises(ValueError, match="descent"):
+ opt._backtracking_line_search(quadratic_objective, x, quadratic_objective.g_x(x))
+
+
+@pytest.mark.parametrize(
+ "method, iterations",
+ [
+ (opt.BFGS, 300),
+ (opt.quasi_newton, 300),
+ (opt.fletcher_reeves, 500),
+ ],
+)
+def test_first_order_methods_reduce_rosenbrock(method, iterations):
+ objective = func.rosenbrock()
+ start = np.array([-1.5, -4.0])
+ path = method(start, objective, iteration=iterations)
+ assert objective.f_x(path[:, -1]) < 1e-8
+ assert np.linalg.norm(objective.g_x(path[:, -1])) < 1e-3
+
+
+def test_figure_generation_preserves_original_outputs(tmp_path):
+ main.generate_figures(tmp_path)
+ expected_sizes = {
+ "quadratic.png": (800, 600),
+ "quasi_newton.png": (800, 600),
+ "cg.png": (1200, 1020),
+ }
+ for filename, expected_size in expected_sizes.items():
+ output = Path(tmp_path, filename)
+ assert output.is_file()
+ assert output.stat().st_size > 10_000
+ with Image.open(output) as image:
+ assert image.size == expected_size