Improve bounds processing if no bounds are specified (#618)

Author: timmens

src/optimagic/differentiation/generate_steps.py

@@ -3,6 +3,8 @@
 
 import numpy as np
 
+from optimagic.utilities import fast_numpy_full
+
 
 class Steps(NamedTuple):
     pos: np.ndarray
@@ -90,12 +92,21 @@ def generate_steps(
     )
     min_steps = base_steps if min_steps is None else min_steps
 
-    assert (bounds.upper - bounds.lower >= 2 * min_steps).all(), (
+    lower_bounds = bounds.lower
+    upper_bounds = bounds.upper
+    # None-valued bounds are handled by instantiating them as an -inf and inf array. In
+    # the future, this should be handled more gracefully.
+    if lower_bounds is None:
+        lower_bounds = fast_numpy_full(len(x), fill_value=-np.inf)
+    if upper_bounds is None:
+        upper_bounds = fast_numpy_full(len(x), fill_value=np.inf)
+
+    assert (upper_bounds - lower_bounds >= 2 * min_steps).all(), (
         "min_steps is too large to fit into bounds."
     )
 
-    upper_step_bounds = bounds.upper - x
-    lower_step_bounds = bounds.lower - x
+    upper_step_bounds = upper_bounds - x
+    lower_step_bounds = lower_bounds - x
 
     pos = step_ratio ** np.arange(n_steps) * base_steps.reshape(-1, 1)
     neg = -pos.copy()
@@ -105,7 +116,7 @@ def generate_steps(
             x, pos, neg, method, lower_step_bounds, upper_step_bounds
         )
 
-    if np.isfinite(bounds.lower).any() or np.isfinite(bounds.upper).any():
+    if np.isfinite(lower_bounds).any() or np.isfinite(upper_bounds).any():
         pos, neg = _rescale_to_accomodate_bounds(
             base_steps, pos, neg, lower_step_bounds, upper_step_bounds, min_steps
         )

src/optimagic/optimization/optimize.py

@@ -561,38 +561,33 @@ def _optimize(problem: OptimizationProblem) -> OptimizeResult:
     # Strict checking if bounds are required and infinite values in bounds
     # ==================================================================================
     if problem.algorithm.algo_info.supports_bounds:
-        # we need all_bounds_infinite and some_bounds_infinite to check
-        # for arrays of np.inf and np.-inf until bounds
-        # processing is improved
-        all_bounds_infinite = (
-            np.isinf(internal_params.lower_bounds).all()
-            and np.isinf(internal_params.upper_bounds).all()
-        )
-        some_bounds_infinite = (
-            np.isinf(internal_params.lower_bounds).all()
-            or np.isinf(internal_params.upper_bounds).all()
-        )
-        infinite_values_in_bounds = (
-            np.isinf(internal_params.lower_bounds).any()
-            or np.isinf(internal_params.upper_bounds).any()
-        )
         bounds_missing = (
             internal_params.lower_bounds is None or internal_params.upper_bounds is None
         )
 
-        if problem.algorithm.algo_info.needs_bounds:
-            if bounds_missing or some_bounds_infinite:
-                raise IncompleteBoundsError(
-                    f"Algorithm {problem.algorithm.name} needs finite bounds "
-                    f"for all parameters. "
-                )
-        if not problem.algorithm.algo_info.supports_infinite_bounds:
-            if not all_bounds_infinite:
-                if infinite_values_in_bounds:
-                    raise IncompleteBoundsError(
-                        f"Algorithm {problem.algorithm.name} does not support "
-                        f"infinite values in bounds for parameters. "
-                    )
+        # Check for infinite values in bounds arrays (only possible in mixed cases now)
+        infinite_values_in_bounds = False
+        if internal_params.lower_bounds is not None:
+            infinite_values_in_bounds |= np.isinf(internal_params.lower_bounds).any()
+        if internal_params.upper_bounds is not None:
+            infinite_values_in_bounds |= np.isinf(internal_params.upper_bounds).any()
+
+        # Case 1: Algorithm needs bounds but none provided
+        if problem.algorithm.algo_info.needs_bounds and bounds_missing:
+            raise IncompleteBoundsError(
+                f"Algorithm {problem.algorithm.name} requires bounds for all "
+                "parameters. Please provide finite lower and upper bounds."
+            )
+
+        # Case 2: Algorithm doesn't support infinite bounds but they are present
+        if (
+            not problem.algorithm.algo_info.supports_infinite_bounds
+            and infinite_values_in_bounds
+        ):
+            raise IncompleteBoundsError(
+                f"Algorithm {problem.algorithm.name} does not support infinite bounds. "
+                "Please provide finite bounds for all parameters."
+            )
 
     # ==================================================================================
     # Do some things that require internal parameters or bounds

src/optimagic/optimizers/fides.py

@@ -183,6 +183,9 @@ def fides_internal(
 
     hessian_instance = _create_hessian_updater_from_user_input(hessian_update_strategy)
 
+    lower_bounds = np.full(len(x), -np.inf) if lower_bounds is None else lower_bounds
+    upper_bounds = np.full(len(x), np.inf) if upper_bounds is None else upper_bounds
+
     opt = Optimizer(
         fun=fun_and_jac,
         lb=lower_bounds,

src/optimagic/optimizers/scipy_optimizers.py

@@ -472,12 +472,15 @@ def _solve_internal_problem(
         else:
             tr_solver_options = self.tr_solver_options
 
+        lower_bounds = -np.inf if problem.bounds.lower is None else problem.bounds.lower
+        upper_bounds = np.inf if problem.bounds.upper is None else problem.bounds.upper
+
         raw_res = scipy.optimize.least_squares(
             fun=problem.fun,
             x0=x0,
             # This optimizer does not work with fun_and_jac
             jac=problem.jac,
-            bounds=(problem.bounds.lower, problem.bounds.upper),
+            bounds=(lower_bounds, upper_bounds),
             method="trf",
             max_nfev=self.stopping_maxfun,
             ftol=self.convergence_ftol_rel,
@@ -522,12 +525,15 @@ def _solve_internal_problem(
         else:
             tr_solver_options = self.tr_solver_options
 
+        lower_bounds = -np.inf if problem.bounds.lower is None else problem.bounds.lower
+        upper_bounds = np.inf if problem.bounds.upper is None else problem.bounds.upper
+
         raw_res = scipy.optimize.least_squares(
             fun=problem.fun,
             x0=x0,
             # This optimizer does not work with fun_and_jac
             jac=problem.jac,
-            bounds=(problem.bounds.lower, problem.bounds.upper),
+            bounds=(lower_bounds, upper_bounds),
             method="dogbox",
             max_nfev=self.stopping_maxfun,
             ftol=self.convergence_ftol_rel,
@@ -1204,8 +1210,13 @@ def _get_workers(n_cores, batch_evaluator):
     return out
 
 
-def _get_scipy_bounds(bounds: InternalBounds) -> ScipyBounds:
-    return ScipyBounds(lb=bounds.lower, ub=bounds.upper)
+def _get_scipy_bounds(bounds: InternalBounds) -> ScipyBounds | None:
+    if bounds.lower is None and bounds.upper is None:
+        return None
+
+    lower = bounds.lower if bounds.lower is not None else -np.inf
+    upper = bounds.upper if bounds.upper is not None else np.inf
+    return ScipyBounds(lb=lower, ub=upper)
 
 
 def process_scipy_result_old(scipy_results_obj):

src/optimagic/optimizers/tao_optimizers.py

@@ -152,10 +152,15 @@ def func_tao(tao, x, resid_out):  # noqa: ARG001
         raise ValueError("The initial trust region radius must be > 0.")
     tao.setInitialTrustRegionRadius(trustregion_initial_radius)
 
-    # Add bounds.
-    lower_bounds = _initialise_petsc_array(lower_bounds)
-    upper_bounds = _initialise_petsc_array(upper_bounds)
-    tao.setVariableBounds(lower_bounds, upper_bounds)
+    # Add bounds if provided.
+    if lower_bounds is not None or upper_bounds is not None:
+        if lower_bounds is None:
+            lower_bounds = np.full(len(x), -np.inf)
+        if upper_bounds is None:
+            upper_bounds = np.full(len(x), np.inf)
+        lower_bounds = _initialise_petsc_array(lower_bounds)
+        upper_bounds = _initialise_petsc_array(upper_bounds)
+        tao.setVariableBounds(lower_bounds, upper_bounds)
 
     # Put the starting values into the container and pass them to the optimizer.
     tao.setInitial(_x)
@@ -197,7 +202,8 @@ def func_tao(tao, x, resid_out):  # noqa: ARG001
     results = _process_pounders_results(residuals_out, tao)
 
     # Destroy petsc objects for memory reasons.
-    for obj in [tao, _x, residuals_out, lower_bounds, upper_bounds]:
+    petsc_bounds = [b for b in (lower_bounds, upper_bounds) if b is not None]
+    for obj in [tao, _x, residuals_out, *petsc_bounds]:
         obj.destroy()
 
     return results

src/optimagic/parameters/bounds.py

@@ -1,7 +1,7 @@
 from __future__ import annotations
 
 from dataclasses import dataclass
-from typing import Any, Sequence
+from typing import Any, Literal, Sequence
 
 import numpy as np
 from numpy.typing import NDArray
@@ -12,6 +12,7 @@
 from optimagic.exceptions import InvalidBoundsError
 from optimagic.parameters.tree_registry import get_registry
 from optimagic.typing import PyTree, PyTreeRegistry
+from optimagic.utilities import fast_numpy_full
 
 
 @dataclass(frozen=True)
@@ -60,8 +61,8 @@ def pre_process_bounds(
 
 
 def _process_bounds_sequence(bounds: Sequence[tuple[float, float]]) -> Bounds:
-    lower = np.full(len(bounds), -np.inf)
-    upper = np.full(len(bounds), np.inf)
+    lower = fast_numpy_full(len(bounds), fill_value=-np.inf)
+    upper = fast_numpy_full(len(bounds), fill_value=np.inf)
 
     for i, (lb, ub) in enumerate(bounds):
         if lb is not None:
@@ -76,14 +77,14 @@ def get_internal_bounds(
     bounds: Bounds | None = None,
     registry: PyTreeRegistry | None = None,
     add_soft_bounds: bool = False,
-) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
+) -> tuple[NDArray[np.float64] | None, NDArray[np.float64] | None]:
     """Create consolidated and flattened bounds for params.
 
     If params is a DataFrame with value column, the user provided bounds are
     extended with bounds from the params DataFrame.
 
-    If no bounds are available the entry is set to minus np.inf for the lower bound and
-    np.inf for the upper bound.
+    If no bounds are provided, we return None. If some bounds are available the missing
+    entries are set to -np.inf for the lower bound and np.inf for the upper bound.
 
     The bounds provided in `bounds` override bounds provided in params if both are
     specified (in the case where params is a DataFrame with bounds as a column).
@@ -109,10 +110,11 @@ def get_internal_bounds(
         add_soft_bounds=add_soft_bounds,
     )
     if fast_path:
-        return _get_fast_path_bounds(
-            params=params,
-            bounds=bounds,
-        )
+        return _get_fast_path_bounds(bounds)
+
+    # Handling of None-valued bounds in the slow path needs to be improved. Currently,
+    # None-valued bounds are replaced with arrays of np.inf and -np.inf, and then
+    # translated back to None if all entries are non-finite.
 
     registry = get_registry(extended=True) if registry is None else registry
     n_params = len(tree_leaves(params, registry=registry))
@@ -149,11 +151,18 @@ def get_internal_bounds(
         msg = "Invalid bounds. Some lower bounds are larger than upper bounds."
         raise InvalidBoundsError(msg)
 
+    if np.isinf(lower_flat).all():
+        lower_flat = None  # type: ignore[assignment]
+    if np.isinf(upper_flat).all():
+        upper_flat = None  # type: ignore[assignment]
+
     return lower_flat, upper_flat
 
 
 def _update_bounds_and_flatten(
-    nan_tree: PyTree, bounds: PyTree, kind: str
+    nan_tree: PyTree,
+    bounds: PyTree,
+    kind: Literal["lower_bound", "upper_bound", "soft_lower_bound", "soft_upper_bound"],
 ) -> NDArray[np.float64]:
     """Flatten bounds array and update it with bounds from params.
 
@@ -213,7 +222,7 @@ def _is_fast_path(params: PyTree, bounds: Bounds, add_soft_bounds: bool) -> bool
     if not _is_1d_array(params):
         out = False
 
-    for bound in bounds.lower, bounds.upper:
+    for bound in (bounds.lower, bounds.upper):
         if not (_is_1d_array(bound) or bound is None):
             out = False
     return out
@@ -224,21 +233,27 @@ def _is_1d_array(candidate: Any) -> bool:
 
 
 def _get_fast_path_bounds(
-    params: PyTree, bounds: Bounds
-) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
+    bounds: Bounds,
+) -> tuple[NDArray[np.float64] | None, NDArray[np.float64] | None]:
     if bounds.lower is None:
-        # faster than np.full
-        lower_bounds = np.array([-np.inf] * len(params))
+        lower_bounds = None
     else:
         lower_bounds = bounds.lower.astype(float)
+        if np.isinf(lower_bounds).all():
+            lower_bounds = None
 
     if bounds.upper is None:
-        # faster than np.full
-        upper_bounds = np.array([np.inf] * len(params))
+        upper_bounds = None
     else:
         upper_bounds = bounds.upper.astype(float)
-
-    if (lower_bounds > upper_bounds).any():
+        if np.isinf(upper_bounds).all():
+            upper_bounds = None
+
+    if (
+        lower_bounds is not None
+        and upper_bounds is not None
+        and (lower_bounds > upper_bounds).any()
+    ):
         msg = "Invalid bounds. Some lower bounds are larger than upper bounds."
         raise InvalidBoundsError(msg)
 

src/optimagic/parameters/consolidate_constraints.py

@@ -11,7 +11,10 @@
 import pandas as pd
 
 from optimagic.exceptions import InvalidConstraintError
-from optimagic.utilities import number_of_triangular_elements_to_dimension
+from optimagic.utilities import (
+    fast_numpy_full,
+    number_of_triangular_elements_to_dimension,
+)
 
 
 def consolidate_constraints(
@@ -25,8 +28,8 @@ def consolidate_constraints(
             constraints have been rewritten as linear constraints and
             pairwise_equality constraints have been rewritten as equality constraints.
         parvec (np.ndarray): 1d numpy array with parameters.
-        lower_bounds (np.ndarray): 1d numpy array with lower_bounds
-        upper_bounds (np.ndarray): 1d numpy array wtih upper_bounds
+        lower_bounds (np.ndarray | None): 1d numpy array with lower_bounds
+        upper_bounds (np.ndarray | None): 1d numpy array with upper_bounds
         param_names (list): Names of parameters. Used for error messages.
 
     Returns:
@@ -37,6 +40,13 @@ def consolidate_constraints(
             constraints.
 
     """
+    # None-valued bounds are handled by instantiating them as an -inf and inf array. In
+    # the future, this should be handled more gracefully.
+    if lower_bounds is None:
+        lower_bounds = fast_numpy_full(len(parvec), fill_value=-np.inf)
+    if upper_bounds is None:
+        upper_bounds = fast_numpy_full(len(parvec), fill_value=np.inf)
+
     raw_eq, other_constraints = _split_constraints(constraints, "equality")
     equality_constraints = _consolidate_equality_constraints(raw_eq)
 

src/optimagic/parameters/conversion.py

@@ -179,12 +179,12 @@ def _fast_derivative_to_internal(
     )
 
     if bounds is None or bounds.lower is None:
-        lower_bounds = np.full(len(params), -np.inf)
+        lower_bounds = None
     else:
         lower_bounds = bounds.lower.astype(float)
 
     if bounds is None or bounds.upper is None:
-        upper_bounds = np.full(len(params), np.inf)
+        upper_bounds = None
     else:
         upper_bounds = bounds.upper.astype(float)
 

src/optimagic/parameters/process_constraints.py

@@ -42,8 +42,8 @@ def process_constraints(
             have already been consolidated into an ``"index"`` field that selects
             the same parameters from the flattened_parameter vector.
         params_vec (np.ndarray): Flattened version of params.
-        lower_bounds (np.ndarray): Lower bounds for params_vec.
-        upper_bounds (np.ndarray): Upper bounds for params_vec.
+        lower_bounds (np.ndarray | None): Lower bounds for params_vec.
+        upper_bounds (np.ndarray | None): Upper bounds for params_vec.
         param_names (list): Names of the flattened parameters. Only used to produce
             good error messages.