Add Hessian check for 2-parameter problems (#363)

* Add Hessian check

* Update test_classification.py

* Update to implicitly concatenated strings

* Update test_classification

* Add tests on insensitivity

* Update CHANGELOG.md

* Fix typo

Co-authored-by: Brady Planden <[email protected]>

* Update input to OptimisationResult

* Improve insensitivity tests

* Improve correlation checks

* Update in line with OptimisationResult

* Make one test a unit test

* Update function name

* Check bounds proximity if infinite cost

---------

Co-authored-by: Brady Planden <[email protected]>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+- [#362](https://github.com/pybop-team/PyBOP/issues/362) - Adds the `classify_using_Hessian` functionality to classify the optimised result.
 - [#584](https://github.com/pybop-team/PyBOP/pull/584) - Adds the `GroupedSPMe` model for parameter identification.
 - [#571](https://github.com/pybop-team/PyBOP/pull/571) - Adds Multistart functionality to optimisers via initialisation arg `multistart`.
 - [#582](https://github.com/pybop-team/PyBOP/pull/582) - Fixes `population_size` arg for Pints' based optimisers, reshapes `parameters.rvs` to be parameter instances.

pybop/__init__.py

@@ -170,6 +170,11 @@
 from .observers.unscented_kalman import UnscentedKalmanFilterObserver
 from .observers.observer import Observer
 
+#
+# Classification classes
+#
+from ._classification import classify_using_hessian
+
 #
 # Plotting classes
 #

pybop/_classification.py

@@ -0,0 +1,139 @@
+from typing import Optional
+
+import numpy as np
+
+from pybop import OptimisationResult
+
+
+def classify_using_hessian(
+    result: OptimisationResult, dx=None, cost_tolerance: Optional[float] = 1e-5
+):
+    """
+    A simple check for parameter correlations based on numerical approximation
+    of the Hessian matrix at the optimal point using central finite differences.
+
+    Parameters
+    ---------
+    result : OptimisationResult
+        The PyBOP optimisation results.
+    dx : array-like, optional
+        An array of small positive values used to check proximity to the parameter
+        bounds and as the perturbation distance in the finite difference calculations.
+    cost_tolerance : float, optional
+        A small positive tolerance used for cost value comparisons (default: 1e-5).
+    """
+    x = result.x
+    dx = np.asarray(dx) if dx is not None else np.maximum(x, 1e-40) * 1e-2
+    final_cost = result.final_cost
+    cost = result.cost
+    parameters = cost.parameters
+    minimising = result.minimising
+
+    n = len(x)
+    if n != 2 or len(dx) != n:
+        raise ValueError(
+            "The function classify_using_hessian currently only works in the case "
+            "of 2 parameters, and dx must have the same length as x."
+        )
+
+    # Get a list of parameter names for use in the output message
+    names = list(parameters.keys())
+
+    # Evaluate the cost for a grid of surrounding points
+    costs = np.zeros((3, 3, 2))
+    for i in np.arange(0, 3):
+        for j in np.arange(0, 3):
+            if i == j == 1:
+                costs[1, 1, 0] = final_cost
+                costs[1, 1, 1] = final_cost
+            else:
+                costs[i, j, 0] = cost(x + np.multiply([i - 1, j - 1], dx))
+                costs[i, j, 1] = cost(x + np.multiply([i - 1, j - 1], 2 * dx))
+
+    def check_proximity_to_bounds(parameters, x, dx, names) -> str:
+        bounds = parameters.get_bounds()
+        if bounds is not None:
+            for i, value in enumerate(x):
+                if value > bounds["upper"][i] - dx[i]:
+                    return f" The result is near the upper bound of {names[i]}."
+
+                if value < bounds["lower"][i] + dx[i]:
+                    return f" The result is near the lower bound of {names[i]}."
+        return ""
+
+    # Classify the result
+    if (minimising and np.any(costs < final_cost)) or (
+        not minimising and np.any(costs > final_cost)
+    ):
+        message = "The optimiser has not converged to a stationary point."
+        message += check_proximity_to_bounds(parameters, x, dx, names)
+
+    elif not np.all([np.isfinite(cost) for cost in costs]):
+        message = "Classification cannot proceed due to infinite cost value(s)."
+        message += check_proximity_to_bounds(parameters, x, dx, names)
+
+    else:
+        # Estimate the Hessian using second-order accurate central finite differences
+        # cfd_hessian = np.zeros((2, 2))
+        # cfd_hessian[0, 0] = costs[2,1,0] - 2 * costs[1,1,0] + costs[0,1,0]
+        # cfd_hessian[0, 1] = (costs[2,2,0] - costs[2,0,0] + costs[0,0,0] - costs[0,2,0]) / 4
+        # cfd_hessian[1, 0] = cfd_hessian[0, 1]
+        # cfd_hessian[1, 1] = costs[1,2,0] - 2 * costs[1,1,0] + costs[1,0,0]
+
+        # Estimate the Hessian using fourth-order accurate central finite differences
+        cfd_hessian = np.zeros((2, 2))
+        cfd_hessian[0, 0] = (
+            -costs[2, 1, 1]
+            + 16 * costs[2, 1, 0]
+            - 30 * costs[1, 1, 0]
+            + 16 * costs[0, 1, 0]
+            - costs[0, 1, 1]
+        ) / 12
+        cfd_hessian[0, 1] = (
+            -(costs[2, 2, 1] - costs[2, 0, 1] + costs[0, 0, 1] - costs[0, 2, 1])
+            + 16 * (costs[2, 2, 0] - costs[2, 0, 0] + costs[0, 0, 0] - costs[0, 2, 0])
+        ) / 48
+        cfd_hessian[1, 0] = cfd_hessian[0, 1]
+        cfd_hessian[1, 1] = (
+            -costs[1, 2, 1]
+            + 16 * costs[1, 2, 0]
+            - 30 * costs[1, 1, 0]
+            + 16 * costs[1, 0, 0]
+            - costs[1, 0, 1]
+        ) / 12
+
+        # Compute the eigenvalues and sort into ascending order
+        eigenvalues, eigenvectors = np.linalg.eig(cfd_hessian)
+        idx = eigenvalues.argsort()
+        eigenvalues = eigenvalues[idx]
+        eigenvectors = eigenvectors[:, idx]
+
+        # Classify the result
+        if np.all(eigenvalues > cost_tolerance):
+            message = "The optimiser has located a minimum."
+        elif np.all(eigenvalues < -cost_tolerance):
+            message = "The optimiser has located a maximum."
+        elif np.all(np.abs(eigenvalues) > cost_tolerance):
+            message = "The optimiser has located a saddle point."
+        elif np.all(np.abs(eigenvalues) < cost_tolerance):
+            message = f"The cost variation is smaller than the cost tolerance: {cost_tolerance}."
+        else:
+            # One eigenvalue is too small to classify with certainty
+            message = "The cost variation is too small to classify with certainty."
+
+        # Check for parameter correlations
+        if np.any(np.abs(eigenvalues) > cost_tolerance):
+            if np.allclose(eigenvectors[0], np.array([1, 0])):
+                message += f" The cost is insensitive to a change of {dx[0]:.2g} in {names[0]}."
+            elif np.allclose(eigenvectors[0], np.array([0, 1])):
+                message += f" The cost is insensitive to a change of {dx[1]:.2g} in {names[1]}."
+            else:
+                diagonal_costs = [
+                    cost(x - np.multiply(eigenvectors[:, 0], dx)),
+                    cost(x + np.multiply(eigenvectors[:, 0], dx)),
+                ]
+                if np.allclose(final_cost, diagonal_costs, atol=cost_tolerance, rtol=0):
+                    message += " There may be a correlation between these parameters."
+
+    print(message)
+    return message

tests/integration/test_classification.py

@@ -0,0 +1,174 @@
+import numpy as np
+import pytest
+from pybamm import Parameter
+
+import pybop
+
+
+class TestClassification:
+    """
+    A class to test the classification of different optimisation results.
+    """
+
+    @pytest.fixture(
+        params=[
+            np.asarray([0.05, 0.05]),
+            np.asarray([0.1, 0.05]),
+            np.asarray([0.05, 0.01]),
+        ]
+    )
+    def parameters(self, request):
+        ground_truth = request.param
+        return pybop.Parameters(
+            pybop.Parameter(
+                "R0 [Ohm]",
+                prior=pybop.Gaussian(0.05, 0.01),
+                bounds=[0.02, 0.08],
+                true_value=ground_truth[0],
+            ),
+            pybop.Parameter(
+                "R1 [Ohm]",
+                prior=pybop.Gaussian(0.05, 0.01),
+                bounds=[0.02, 0.08],
+                true_value=ground_truth[1],
+            ),
+        )
+
+    @pytest.fixture
+    def parameter_set(self):
+        parameter_set = pybop.ParameterSet(
+            json_path="examples/parameters/initial_ecm_parameters.json"
+        )
+        parameter_set.import_parameters()
+        parameter_set.params.update({"C1 [F]": 1000})
+        return parameter_set
+
+    @pytest.fixture
+    def model(self, parameter_set, parameters):
+        parameter_set.params.update(parameters.as_dict(parameters.true_value()))
+        return pybop.empirical.Thevenin(parameter_set=parameter_set)
+
+    @pytest.fixture
+    def dataset(self, model):
+        experiment = pybop.Experiment(
+            [
+                "Discharge at 0.5C for 2 minutes (4 seconds period)",
+                "Charge at 0.5C for 2 minutes (4 seconds period)",
+            ]
+        )
+        solution = model.predict(experiment=experiment)
+        return pybop.Dataset(
+            {
+                "Time [s]": solution["Time [s]"].data,
+                "Current function [A]": solution["Current [A]"].data,
+                "Voltage [V]": solution["Voltage [V]"].data,
+            }
+        )
+
+    @pytest.fixture
+    def problem(self, model, parameters, dataset):
+        return pybop.FittingProblem(model, parameters, dataset)
+
+    @pytest.mark.integration
+    def test_classify_using_hessian(self, problem):
+        cost = pybop.RootMeanSquaredError(problem)
+        x = cost.parameters.true_value()
+        bounds = cost.parameters.get_bounds()
+        x0 = np.clip(x, bounds["lower"], bounds["upper"])
+        optim = pybop.Optimisation(cost=cost)
+        results = pybop.OptimisationResult(x=x0, optim=optim)
+
+        if np.all(x == np.asarray([0.05, 0.05])):
+            message = pybop.classify_using_hessian(results)
+            assert message == "The optimiser has located a minimum."
+        elif np.all(x == np.asarray([0.1, 0.05])):
+            message = pybop.classify_using_hessian(results)
+            assert message == (
+                "The optimiser has not converged to a stationary point."
+                " The result is near the upper bound of R0 [Ohm]."
+            )
+        elif np.all(x == np.asarray([0.05, 0.01])):
+            message = pybop.classify_using_hessian(results)
+            assert message == (
+                "The optimiser has not converged to a stationary point."
+                " The result is near the lower bound of R1 [Ohm]."
+            )
+        else:
+            raise Exception(f"Please add a check for these values: {x}")
+
+        if np.all(x == np.asarray([0.05, 0.05])):
+            cost = pybop.GaussianLogLikelihoodKnownSigma(problem, sigma0=0.002)
+            optim = pybop.Optimisation(cost=cost)
+            results = pybop.OptimisationResult(x=x, optim=optim)
+
+            message = pybop.classify_using_hessian(results)
+            assert message == "The optimiser has located a maximum."
+
+        # message = pybop.classify_using_hessian(results)
+        # assert message == "The optimiser has located a saddle point."
+
+    @pytest.mark.integration
+    def test_insensitive_classify_using_hessian(self, parameter_set):
+        param_R0_a = pybop.Parameter(
+            "R0_a [Ohm]",
+            bounds=[0, 0.002],
+            true_value=0.001,
+        )
+        param_R0_b = pybop.Parameter(
+            "R0_b [Ohm]",
+            bounds=[-1e-4, 1e-4],
+            true_value=0,
+        )
+        parameter_set.params.update(
+            {"R0_a [Ohm]": 0.001, "R0_b [Ohm]": 0},
+            check_already_exists=False,
+        )
+        parameter_set.params.update(
+            {"R0 [Ohm]": Parameter("R0_a [Ohm]") + Parameter("R0_b [Ohm]")},
+        )
+        model = pybop.empirical.Thevenin(parameter_set=parameter_set)
+
+        experiment = pybop.Experiment(
+            ["Discharge at 0.5C for 2 minutes (4 seconds period)"]
+        )
+        solution = model.predict(experiment=experiment)
+        dataset = pybop.Dataset(
+            {
+                "Time [s]": solution["Time [s]"].data,
+                "Current function [A]": solution["Current [A]"].data,
+                "Voltage [V]": solution["Voltage [V]"].data,
+            }
+        )
+
+        for parameters in [
+            pybop.Parameters(param_R0_b, param_R0_a),
+            pybop.Parameters(param_R0_a, param_R0_b),
+        ]:
+            problem = pybop.FittingProblem(model, parameters, dataset)
+            cost = pybop.SumofPower(problem, p=1)
+            x = cost.parameters.true_value()
+            optim = pybop.Optimisation(cost=cost)
+            results = pybop.OptimisationResult(x=x, optim=optim)
+
+            message = pybop.classify_using_hessian(results)
+            assert message == (
+                "The cost variation is too small to classify with certainty."
+                " The cost is insensitive to a change of 1e-42 in R0_b [Ohm]."
+            )
+
+        message = pybop.classify_using_hessian(results, dx=[0.0001, 0.0001])
+        assert message == (
+            "The optimiser has located a minimum."
+            " There may be a correlation between these parameters."
+        )
+
+        message = pybop.classify_using_hessian(results, cost_tolerance=1e-2)
+        assert message == (
+            "The cost variation is smaller than the cost tolerance: 0.01."
+        )
+
+        message = pybop.classify_using_hessian(results, dx=[1, 1])
+        assert message == (
+            "Classification cannot proceed due to infinite cost value(s)."
+            " The result is near the upper bound of R0_a [Ohm]."
+        )

tests/unit/test_classifier.py

@@ -0,0 +1,50 @@
+import numpy as np
+import pytest
+
+import pybop
+
+
+class TestClassifier:
+    """
+    A class to test the classification of different optimisation results.
+    """
+
+    @pytest.fixture
+    def problem(self):
+        model = pybop.empirical.Thevenin()
+        experiment = pybop.Experiment(
+            [
+                "Discharge at 0.5C for 2 minutes (4 seconds period)",
+                "Charge at 0.5C for 2 minutes (4 seconds period)",
+            ]
+        )
+        solution = model.predict(experiment=experiment)
+        dataset = pybop.Dataset(
+            {
+                "Time [s]": solution["Time [s]"].data,
+                "Current function [A]": solution["Current [A]"].data,
+                "Voltage [V]": solution["Voltage [V]"].data,
+            }
+        )
+        parameters = pybop.Parameters(
+            pybop.Parameter(
+                "R0 [Ohm]",
+                prior=pybop.Uniform(0.001, 0.1),
+                bounds=[1e-4, 0.1],
+            ),
+        )
+        return pybop.FittingProblem(model, parameters, dataset)
+
+    @pytest.mark.unit
+    def test_classify_using_hessian_invalid(self, problem):
+        cost = pybop.SumSquaredError(problem)
+        optim = pybop.Optimisation(cost=cost)
+        x = np.asarray([0.001])
+        results = pybop.OptimisationResult(x=x, optim=optim)
+
+        with pytest.raises(
+            ValueError,
+            match="The function classify_using_hessian currently only works"
+            " in the case of 2 parameters, and dx must have the same length as x.",
+        ):
+            pybop.classify_using_hessian(results)