Repository: OpenSourceEconomics/estimagic

Issue:
Allow user to specify steps for numdifftools, scaling
I set the Richardson extrapolation version running yesterday. Still running.

I would be willing to sacrifice some precision for speed here; surely for testing. So it would be great if it was possible to pass the options `step_ratio=2.0, step=1, order=1, num_terms=2` there. Although I must admit that I have not gotten my head what precisely they do yet.

I have not really verified, but it seems to me that the steps are absolute. They certainly were for the manually specified finite-differencing derivative. With large parameter values, this leads to numerical error dominating. To be on the safe side here, I would suggest dividing by the absolute value of the parameter value.

Discussion:

Patch:
diff --git a/docs/source/differentiation.rst b/docs/source/differentiation.rst
index fb7a17d..74d0fa0 100644
--- a/docs/source/differentiation.rst
+++ b/docs/source/differentiation.rst
@@ -15,8 +15,6 @@ backward or central differences for the following problems:
 - Jacobian matrix of real vector-valued multivariate functions
 - Hessian matrix of real valued multivariate functions
 
-The methods explained below are not implemented by us and instead called from the
-
 
 Functions
 ---------
@@ -31,9 +29,9 @@ The variable method "method" allows the user to call a specific finite differenc
 method. We only allow for one of the methods and therefore exclude complex step and
 other possibilities:
 
-- forward differences ("method" = "forward")
-- backward differences ("method" = "backward")
-- central differences ("method" = "central")
+- forward differences (method="forward")
+- backward differences (method="backward")
+- central differences (method="central")
 
 As the central difference method is the most accurate, it is the default for all three
 functions. The numdifftools functions use Richardson extrapolation as their standard
@@ -60,6 +58,43 @@ of the jacobian matrix:
 .. autofunction:: jacobian
 
 
+.. _step_options:
+
+Influencing the Step Size
+-------------------------
+
+As mentioned before, we use numdifftools_ for all numerical differentiations with
+Richardson extrapolations. numdifftools_ offers many ways of influencing the step
+size of the extrapolation:
+
+
+- base_step (float, array-like, optional):
+    Defines the minimum step, if None, the value is set to ``EPS**(1/scale)````
+
+- step_ratio (float , optional):
+    Ratio between sequential steps generated. Must be > 1. If None
+    then ``step_ratio`` is 2 for first order derivatives, otherwise it is 1.6
+
+- num_steps (scalar integer, optional):
+    default ``min_num_steps + num_extrap``. Defines number of steps generated.
+    It should be larger than ``min_num_steps = (n + order - 1) / fact`` where
+    ``fact`` is 1, 2 or 4 depending on differentiation method used.
+
+
+- offset (float, optional):
+    offset to the base step.
+
+- num_extrap (int):
+    number of points used for extrapolation.  Numdifftools Documentation says
+    that the default is 0 but that would be surprising.
+
+
+- scale (float, array like):
+    scale used in base step. If not None it will override the default
+    computed with the default_scale function.
+
+
+
 Methods
 -------
 
diff --git a/environment.yml b/environment.yml
index 77da094..e671178 100644
--- a/environment.yml
+++ b/environment.yml
@@ -32,3 +32,4 @@ dependencies:
     - restructuredtext_lint
     - tox-conda
     - black
+    - pdbpp
diff --git a/estimagic/differentiation/differentiation.py b/estimagic/differentiation/differentiation.py
index 5679cf8..080e397 100644
--- a/estimagic/differentiation/differentiation.py
+++ b/estimagic/differentiation/differentiation.py
@@ -6,7 +6,13 @@ from estimagic.differentiation import differentiation_auxiliary as aux
 
 
 def gradient(
-    func, params, method="central", extrapolation=True, func_args=None, func_kwargs=None
+    func,
+    params,
+    method="central",
+    extrapolation=True,
+    func_args=None,
+    func_kwargs=None,
+    step_options=None,
 ):
     """
     Calculate the gradient of *func*.
@@ -20,11 +26,16 @@ def gradient(
             richardson extrapolation.
         func_args (list): additional positional arguments for func.
         func_kwargs (dict): additional positional arguments for func.
+        step_options (dict): Options for the numdifftools step generator.
+            See :ref:`step_options`
+
 
     Returns:
         Series: The index is the index of params_sr.
 
     """
+    step_options = step_options if step_options is not None else {}
+
     if method not in ["central", "forward", "backward"]:
         raise ValueError("Method has to be in ['central', 'forward', 'backward']")
 
@@ -35,7 +46,9 @@ def gradient(
     params_value = params["value"].to_numpy()
 
     if extrapolation:
-        grad_np = nd.Gradient(internal_func, method=method)(params_value)
+        grad_np = nd.Gradient(internal_func, method=method, **step_options)(
+            params_value
+        )
     else:
         grad_np = _no_extrapolation_gradient(internal_func, params_value, method)
     return pd.Series(data=grad_np, index=params.index, name="gradient")
@@ -52,7 +65,13 @@ def _no_extrapolation_gradient(internal_func, params_value, method):
 
 
 def jacobian(
-    func, params, method="central", extrapolation=True, func_args=None, func_kwargs=None
+    func,
+    params,
+    method="central",
+    extrapolation=True,
+    func_args=None,
+    func_kwargs=None,
+    step_options=None,
 ):
     """
     Calculate the jacobian of *func*.
@@ -68,6 +87,8 @@ def jacobian(
                                 richardson extrapolation.
         func_args (list): additional positional arguments for func.
         func_kwargs (dict): additional positional arguments for func.
+        step_options (dict): Options for the numdifftools step generator.
+            See :ref:`step_options`
 
     Returns:
         DataFrame: If func returns a Series, the index is the index of this Series or
@@ -75,6 +96,8 @@ def jacobian(
         index of params_sr.
 
     """
+    step_options = step_options if step_options is not None else {}
+
     if method not in ["central", "forward", "backward"]:
         raise ValueError("Method has to be in ['central', 'forward', 'backward']")
 
@@ -87,7 +110,7 @@ def jacobian(
     params_value = params["value"].to_numpy()
 
     if extrapolation:
-        jac_np = nd.Jacobian(internal_func, method=method)(params_value)
+        jac_np = nd.Jacobian(internal_func, method=method, **step_options)(params_value)
     else:
         jac_np = _no_extrapolation_jacobian(internal_func, params_value, method)
 
@@ -110,7 +133,13 @@ def _no_extrapolation_jacobian(internal_func, params_value, method):
 
 
 def hessian(
-    func, params, method="central", extrapolation=True, func_args=None, func_kwargs=None
+    func,
+    params,
+    method="central",
+    extrapolation=True,
+    func_args=None,
+    func_kwargs=None,
+    step_options=None,
 ):
     """
     Calculate the hessian of *func*.
@@ -124,11 +153,15 @@ def hessian(
                                 richardson extrapolation.
         func_args (list): additional positional arguments for func.
         func_kwargs (dict): additional positional arguments for func.
+        step_options (dict): Options for the numdifftools step generator.
+            See :ref:`step_options`
 
     Returns:
         DataFrame: The index and columns are the index of params_sr.
 
     """
+    step_options = step_options if step_options is not None else {}
+
     if method != "central":
         raise ValueError("Only the method 'central' is supported.")
 
@@ -139,7 +172,7 @@ def hessian(
     params_value = params["value"].to_numpy()
 
     if extrapolation:
-        hess_np = nd.Hessian(internal_func, method=method)(params_value)
+        hess_np = nd.Hessian(internal_func, method=method, **step_options)(params_value)
     else:
         hess_np = _no_extrapolation_hessian(internal_func, params_value, method)
     return pd.DataFrame(data=hess_np, index=params.index, columns=params.index)
diff --git a/estimagic/optimization/optimize.py b/estimagic/optimization/optimize.py
index 98dabd7..b08f35e 100644
--- a/estimagic/optimization/optimize.py
+++ b/estimagic/optimization/optimize.py
@@ -177,7 +177,7 @@ def minimize(
         )
         outer_server_process.start()
 
-    result, timing_info = _minimize(
+    result = _minimize(
         criterion=criterion,
         criterion_args=criterion_args,
         criterion_kwargs=criterion_kwargs,
@@ -194,7 +194,7 @@ def minimize(
     if dashboard:
         stop_signal.set()
         outer_server_process.terminate()
-    return result, timing_info
+    return result
 
 
 def _minimize(
diff --git a/setup.py b/setup.py
index 52a2736..2b52d31 100644
--- a/setup.py
+++ b/setup.py
@@ -1,7 +1,7 @@
 from setuptools import find_packages
 from setuptools import setup
 
-version = "0.0.15"
+version = "0.0.16"
 
 setup(
     name="estimagic",


Test patch:
diff --git a/estimagic/tests/differentiation/test_differentiation.py b/estimagic/tests/differentiation/test_differentiation.py
index c114685..47c5d54 100644
--- a/estimagic/tests/differentiation/test_differentiation.py
+++ b/estimagic/tests/differentiation/test_differentiation.py
@@ -35,8 +35,9 @@ def test_gradient(statsmodels_fixtures, method, extrapolation):
         logit_loglike,
         fix["params"],
         method="forward",
-        extrapolation=False,
+        extrapolation=extrapolation,
         func_args=func_args,
+        step_options={"step_ratio": 2.0},
     )
     expected = fix["gradient"]
     assert_series_equal(calculated, expected)
@@ -52,6 +53,7 @@ def test_jacobian(statsmodels_fixtures, method, extrapolation):
         method=method,
         extrapolation=extrapolation,
         func_kwargs=func_kwargs,
+        step_options={"step_ratio": 2.0},
     )
 
     expected = fix["jacobian"]
@@ -70,6 +72,7 @@ def test_hessian(statsmodels_fixtures, method, extrapolation):
         method=method,
         extrapolation=extrapolation,
         func_kwargs={"y": fix["y"], "x": fix["x"]},
+        step_options={"step_ratio": 2.0},
     )
     expected = fix["hessian"]
     assert_frame_equal(calculated, expected)