Infinite Gradient Handling

src/optimagic/optimization/internal_optimization_problem.py

@@ -11,7 +11,11 @@
 from optimagic.batch_evaluators import process_batch_evaluator
 from optimagic.differentiation.derivatives import first_derivative
 from optimagic.differentiation.numdiff_options import NumdiffOptions
-from optimagic.exceptions import UserFunctionRuntimeError, get_traceback
+from optimagic.exceptions import (
+    InvalidFunctionError,
+    UserFunctionRuntimeError,
+    get_traceback,
+)
 from optimagic.logging.logger import LogStore
 from optimagic.logging.types import IterationState
 from optimagic.optimization.fun_value import (
@@ -471,6 +475,7 @@ def _pure_evaluate_jac(
         out_jac = _process_jac_value(
             value=jac_value, direction=self._direction, converter=self._converter, x=x
         )
+        self._assert_finite_jac(out_jac, jac_value, params)
 
         stop_time = time.perf_counter()
 
@@ -508,6 +513,7 @@ def func(x: NDArray[np.float64]) -> SpecificFunctionValue:
             p = self._converter.params_from_internal(x)
             return self._fun(p)
 
+        params = self._converter.params_from_internal(x)
         try:
             numdiff_res = first_derivative(
                 func,
@@ -543,6 +549,8 @@ def func(x: NDArray[np.float64]) -> SpecificFunctionValue:
                 warnings.warn(msg)
                 fun_value, jac_value = self._error_penalty_func(x)
 
+        self._assert_finite_jac(jac_value, jac_value, params)
+
         algo_fun_value, hist_fun_value = _process_fun_value(
             value=fun_value,  # type: ignore
             solver_type=self._solver_type,
@@ -682,6 +690,8 @@ def _pure_evaluate_fun_and_jac(
         if self._direction == Direction.MAXIMIZE:
             out_jac = -out_jac
 
+        self._assert_finite_jac(out_jac, jac_value, params)
+
         stop_time = time.perf_counter()
 
         hist_entry = HistoryEntry(
@@ -704,6 +714,30 @@ def _pure_evaluate_fun_and_jac(
 
         return (algo_fun_value, out_jac), hist_entry, log_entry
 
+    def _assert_finite_jac(
+        self, out_jac: NDArray[np.float64], jac_value: PyTree, params: PyTree
+    ) -> None:
+        """Check for infinite and NaN values in the jacobian and raise an error if
+        found.
+
+        Args:
+            out_jac: internal processed gradient to check for infinities.
+            jac_value: original gradient value as returned by the user function,
+                    included in error messages for debugging.
+            params: user-facing parameter representation at evaluation point.
+
+        Raises:
+            InvalidFunctionError: If any infinite values are found in the gradient.
+
+        """
+        if not np.all(np.isfinite(out_jac)):
+            msg = (
+                "Infinite or NaN values found in gradient.\n"
+                f"Parameters: {params},\n"
+                f"Gradient: {jac_value}"
+            )
+            raise InvalidFunctionError(msg)
+
 
 def _process_fun_value(
     value: SpecificFunctionValue,

src/optimagic/parameters/space_conversion.py

@@ -145,7 +145,6 @@ def get_space_converter(
         soft_lower_bounds=_soft_lower,
         soft_upper_bounds=_soft_upper,
     )
-
     return converter, params
 
 

tests/optimagic/optimization/test_invalid_jacobian_value.py

@@ -0,0 +1,119 @@
+import numpy as np
+import pandas as pd
+import pytest
+
+from optimagic.exceptions import InvalidFunctionError
+from optimagic.optimization.fun_value import (
+    LeastSquaresFunctionValue,
+    LikelihoodFunctionValue,
+    ScalarFunctionValue,
+)
+from optimagic.optimization.optimize import minimize
+
+SCALAR_VALUES = [
+    ScalarFunctionValue(5),
+]
+
+LS_VALUES = [
+    LeastSquaresFunctionValue(np.array([1, 2])),
+    LeastSquaresFunctionValue({"a": 1, "b": 2}),
+]
+
+LIKELIHOOD_VALUES = [
+    LikelihoodFunctionValue(np.array([1, 4])),
+    LikelihoodFunctionValue({"a": 1, "b": 4}),
+]
+
+
+def test_with_infinite_jacobian_value_in_lists():
+    def sphere(params):
+        return params @ params
+
+    def sphere_gradient(params):
+        grad = 2 * params
+        grad[(abs(grad) < 1.0) & (abs(grad) > 0.0)] = (
+            np.sign(grad)[(abs(grad) < 1.0) & (abs(grad) > 0.0)] * np.inf
+        )
+        return grad
+
+    with pytest.raises(InvalidFunctionError):
+        minimize(
+            fun=sphere,
+            params=np.arange(10) + 400,
+            algorithm="scipy_lbfgsb",
+            jac=sphere_gradient,
+        )
+
+
+def test_with_infinite_jacobian_value_in_dicts():
+    def sphere(params):
+        return params["a"] ** 2 + params["b"] ** 2 + (params["c"] ** 2).sum()
+
+    def sphere_gradient(params):
+        grad = {
+            "a": 2 * params["a"]
+            if not ((abs(params["a"]) < 1.0) & (abs(params["a"]) > 0.0))
+            else np.sign(params["a"]) * np.inf,
+            "b": 2 * params["b"]
+            if not ((abs(params["b"]) < 1.0) & (abs(params["b"]) > 0.0))
+            else np.sign(params["b"]) * np.inf,
+            "c": 2 * params["c"]
+            if not ((abs(params["c"].sum()) < 1.0) & (abs(params["c"].sum()) > 0.0))
+            else np.sign(params["c"]) * np.inf,
+        }
+        return grad
+
+    with pytest.raises(InvalidFunctionError):
+        minimize(
+            fun=sphere,
+            params={"a": 400, "b": 400, "c": pd.Series([200, 300, 400])},
+            algorithm="scipy_lbfgsb",
+            jac=sphere_gradient,
+        )
+
+
+def test_with_nan_jacobian_value_in_lists():
+    def sphere(params):
+        return params @ params
+
+    def sphere_gradient(params):
+        grad = 2 * params
+        grad[(abs(grad) < 1.0) & (abs(grad) > 0.0)] = (
+            np.sign(grad)[(abs(grad) < 1.0) & (abs(grad) > 0.0)] * np.nan
+        )
+        return grad
+
+    with pytest.raises(InvalidFunctionError):
+        minimize(
+            fun=sphere,
+            params=np.arange(10) + 400,
+            algorithm="scipy_lbfgsb",
+            jac=sphere_gradient,
+        )
+
+
+def test_with_nan_jacobian_value_in_dicts():
+    def sphere(params):
+        return params["a"] ** 2 + params["b"] ** 2 + (params["c"] ** 2).sum()
+
+    def sphere_gradient(params):
+        grad = {
+            "a": 2 * params["a"]
+            if not ((abs(params["a"]) < 1.0) & (abs(params["a"]) > 0.0))
+            else np.sign(params["a"]) * np.nan,
+            "b": 2 * params["b"]
+            if not ((abs(params["b"]) < 1.0) & (abs(params["b"]) > 0.0))
+            else np.sign(params["b"]) * np.nan,
+            "c": 2 * params["c"]
+            if not ((abs(params["c"].sum()) < 1.0) & (abs(params["c"].sum()) > 0.0))
+            else np.sign(params["c"]) * np.nan,
+        }
+        return grad
+
+    with pytest.raises(InvalidFunctionError):
+        minimize(
+            fun=sphere,
+            params={"a": 400, "b": 400, "c": pd.Series([200, 300, 400])},
+            algorithm="scipy_lbfgsb",
+            jac=sphere_gradient,
+        )