Fix optimizer bugs and preserve original figures

.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

.gitignore

@@ -33,3 +33,8 @@ nosetests.xml
 .mr.developer.cfg
 .project
 .pydevproject
+
+# Local development
+.venv/
+__pycache__/
+.pytest_cache/

README.md

@@ -1,8 +1,23 @@
-gradient
-========
-
-optimisation practices
--    Newton, modified Newton, Quasi Newton (rank one with backtracking)
-<img src="https://github.com/wyli/gradient/raw/master/quadratic.png">
-<img src="https://github.com/wyli/gradient/raw/master/quasi_newton.png">
-<img src="https://github.com/wyli/gradient/raw/master/cg.png">
+# 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
+```
+
+![Quadratic comparison](quadratic.png)
+![Newton and quasi-Newton comparison](quasi_newton.png)
+![Rosenbrock method comparison](cg.png)

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

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))