Correct BFGS update with outer product

Author: timothy-nunn

src/pyvmcon/vmcon.py

@@ -504,6 +504,18 @@ def _powells_gamma(gamma: np.ndarray, ksi: np.ndarray, B: np.ndarray) -> np.ndar
     return theta * gamma + (1 - theta) * (B @ ksi)  # eqn 9
 
 
+def _revise_B(current_B: np.ndarray, ksi: np.ndarray, gamma: np.ndarray) -> np.ndarray:
+    """Revises B using a BFGS update.
+
+    Implements Equation 8 of the Crane report.
+    """
+    return (
+        current_B
+        - ((current_B @ np.outer(ksi, ksi) @ current_B.T) / (ksi.T @ current_B @ ksi))
+        + (np.outer(gamma, gamma) / (ksi.T @ gamma))
+    )
+
+
 def calculate_new_B(
     result: Result,
     new_result: Result,
@@ -513,10 +525,7 @@ def calculate_new_B(
     lamda_equality: np.ndarray,
     lamda_inequality: np.ndarray,
 ) -> np.ndarray:
-    """Updates the hessian approximation matrix.
-
-    Uses Equation 8 of the Crane report.
-    """
+    """Updates the hessian approximation matrix."""
     # xi (the symbol name) would be a bit confusing in this context,
     # ksi is how its pronounced in modern greek
     # reshape ksi to be a matrix
@@ -542,11 +551,7 @@ def calculate_new_B(
         logger.warning("All xi (ksi) components are 0")
         ksi[:] = 1e-10
 
-    # eqn 8
-    B_ksi = B @ ksi
-    B += (gamma @ gamma.T) / (ksi.T @ gamma) - ((B_ksi @ ksi.T @ B) / (ksi.T @ B_ksi))
-
-    return B
+    return _revise_B(B, ksi, gamma)
 
 
 def _find_out_of_bounds_vars(higher: np.ndarray, lower: np.ndarray) -> list[str]:

tests/test_vmcon.py

@@ -0,0 +1,65 @@
+"""Test individual units of the PyVMCON implementation."""
+
+from dataclasses import dataclass
+
+import numpy as np
+import pytest
+
+from pyvmcon.vmcon import _revise_B
+
+
+@dataclass
+class BRevisionAsset:
+    """Test asset for testing B matrix revision."""
+
+    B: np.ndarray
+    ksi: np.ndarray
+    eta: np.ndarray
+    expected_new_B: np.ndarray  # noqa: N815
+
+
+@pytest.mark.parametrize(
+    "test_asset",
+    [
+        BRevisionAsset(
+            B=np.identity(2),
+            ksi=np.array([-0.66666666666666663, -0.83333333333333348]),
+            eta=np.array([-1.3425925925925923, -1.7129629629629632]),
+            expected_new_B=np.array(
+                [
+                    [1.3858727457706470, 0.50241291449459347],
+                    [0.50241291449459347, 1.6536252239598812],
+                ]
+            ),
+        ),
+        BRevisionAsset(
+            B=np.array(
+                [
+                    [2.1875467668036239, 1.4714414127452644],
+                    [1.4714414127452644, 2.7501870672148332],
+                ]
+            ),
+            ksi=np.array([-1.2385140125071988e-6, -6.1925700625482853e-7]),
+            eta=np.array([-3.6205427119684330e-6, -3.5255433852299234e-6]),
+            expected_new_B=np.array(
+                [
+                    [2.1875592084316073, 1.4714730232083293],
+                    [1.4714730232083293, 2.7502368318126678],
+                ]
+            ),
+        ),
+    ],
+)
+def test_revise_B(test_asset):
+    """Tests the hessian update implementation.
+
+    Uses data from Example 1 of the NEA (Crane) to ensure PyVMCON agrees with that
+    implementation to at least 14 decimal places.
+    """
+    new_B = _revise_B(test_asset.B, test_asset.ksi, test_asset.eta)
+
+    # check symmetric
+    np.testing.assert_array_almost_equal(new_B, new_B.T, decimal=14)
+
+    # check our revision agrees with NEA version of VMCON
+    np.testing.assert_array_almost_equal(new_B, test_asset.expected_new_B, decimal=14)