Repository: OpenSourceEconomics/estimagic

Issue:
Add scaling of the parameter vector
Scale the parameter vector by the inverse of the start parameters in reparametrize_to_internal and undo the scaling in reparametrize_from_internal.

Discussion:
tobiasraabe: @janos There is an issue documented starting in ``internal_func`` in differentiation. It is because for scaling I have to set up an internal criterion function and pass it to gradient. There it gets wrapped again.
janosg: Whether to take derivatives for standard errors on external or internal params was discussed extensively in #57 and we decided that for now we take the derivatives on on the external params. 

Our derivative functions should work for functions that take numpy arrays as well. If they currently don't we have to rewrite `_create_internal_func`. The double wrapping is not very nice but the speed penalty is probably negligible. 

For the scaling we should probably use faster but less precise derivatives. So `extrapolation=False` and maybe even `method="forward"`. We should also use a larger stepsize. This can currently not be specified but I'll change it in #59. 


tobiasraabe: @janosg I implemented the scaling like you suggested. The linting catches some unused arguments. You can go through the list or I can disable it.
janosg: @peisenha  I requested a review from you because you have implemented precondotioning for an optimizer before. 

Patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 39686e2..28af28e 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -34,13 +34,12 @@ repos:
     hooks:
     -   id: black
         types: [python]
+-   repo: https://github.com/PyCQA/doc8
+    rev: 0.8.1rc2
+    hooks:
+        -   id: doc8
 -   repo: local
     hooks:
-    -   id: doc8
-        entry: doc8
-        language: system
-        types: [rst]
-        name: doc8 - Lint the documentation.
     -   id: restructuredtext-lint
         entry: rst-lint
         language: system
diff --git a/docs/source/estimation.rst b/docs/source/estimation.rst
index abc3e08..f2f4e6b 100644
--- a/docs/source/estimation.rst
+++ b/docs/source/estimation.rst
@@ -184,3 +184,29 @@ The options can be found `here
 - scipy_L-BFGS-B
 - scipy_TNC
 - scipy_SLSQP
+
+
+.. _estimation_general_options:
+
+General options
+---------------
+
+You can pass a dictionary with general options to the
+:func:`~estimagic.optimization.optimize.minimize` or
+:func:`~estimagic.optimization.optimize.maximize` function. The following options are
+available.
+
+**Scaling of the parameter vector**
+
+``estimagic`` is able to automatically scale the parameter vector such that the
+optimizer receives parameters closer to 1. There are two available methods.
+
+1. ``{"scaling": "start_values"}`` divides the parameter vector by the starting
+   values for all starting values not in ``[-1, 1]``.
+
+2. ``{"scaling": "gradient"}`` divides the parameter vector by the inverse of the
+   gradient for each parameter not in ``[-1e-2, 1e-2]``. By default,
+   ``{"scaling_gradient_method": "central", "scaling_gradient_extrapolation": False}``
+   are also set. ``"forward"`` and ``"backward"`` can also be used as a method and
+   extrapolation can be turned on by setting the value to ``True``. Note that without
+   extrapolation, the gradient is computed faster.
diff --git a/docs/source/inference.rst b/docs/source/inference.rst
index 466d0e3..4e2be80 100644
--- a/docs/source/inference.rst
+++ b/docs/source/inference.rst
@@ -4,7 +4,8 @@ Inference
 
 Here we will have standard error functions for the following types of models:
 
-- maximum likelihood (ml) generalized method of moments (gmm) method of simulated moments
+- maximum likelihood (ml) generalized method of moments (gmm) method of simulated
+  moments
 - (msm) indirect inference (ii)
 
 Generalized Method of Moments (GMM)
@@ -17,10 +18,10 @@ the notation and formulations provided in section 13.25 of
 
 Bruce E. Hansen - Econometrics (https://www.ssc.wisc.edu/~bhansen/econometrics)
 
-The distribution of the estimator is shaped by the moment conditions, which are functions
-of the estimated parameter :math:`\beta` to the real numbers. We denote the vector of all
-moment condition values for observation :math:`i` by :math:`\textbf{g}_i`. With this on
-hand the calculation of the covariance matrix can be formalized:
+The distribution of the estimator is shaped by the moment conditions, which are
+functions of the estimated parameter :math:`\beta` to the real numbers. We denote the
+vector of all moment condition values for observation :math:`i` by :math:`\textbf{g}_i`.
+With this on hand the calculation of the covariance matrix can be formalized:
 
 .. math:: V = n^{-1}(Q^{T}WQ)^{-1} (Q^{T}W\Omega WQ) (Q^{T}WQ)^{-1}
 
diff --git a/docs/source/optimize.rst b/docs/source/optimize.rst
index e8cf061..bb1bd2d 100644
--- a/docs/source/optimize.rst
+++ b/docs/source/optimize.rst
@@ -2,7 +2,6 @@
 Optimization tools
 ==================
 
-
 .. _minimize:
 
 minimize function
@@ -25,5 +24,4 @@ Processing constraints
 Enforcing constraints by reparametrizations and bounds
 ------------------------------------------------------
 
-
 .. automodule:: estimagic.optimization.reparametrize
diff --git a/environment.yml b/environment.yml
index 266bba2..77da094 100644
--- a/environment.yml
+++ b/environment.yml
@@ -25,7 +25,6 @@ dependencies:
   - jupyter
   - fuzzywuzzy
   - pip:
-    - doc8
     - pre-commit
     - pytest
     - pytest-mock
diff --git a/estimagic/differentiation/differentiation.py b/estimagic/differentiation/differentiation.py
index f736f5d..5679cf8 100644
--- a/estimagic/differentiation/differentiation.py
+++ b/estimagic/differentiation/differentiation.py
@@ -174,6 +174,7 @@ def _create_internal_func(func, params, func_args, func_kwargs):
     def internal_func(x):
         p = params.copy(deep=True)
         p["value"] = x
+
         func_value = func(p, *func_args, **func_kwargs)
         if isinstance(func_value, (pd.DataFrame, pd.Series)):
             func_value = func_value.to_numpy()
diff --git a/estimagic/optimization/optimize.py b/estimagic/optimization/optimize.py
index 4fe942d..98dabd7 100644
--- a/estimagic/optimization/optimize.py
+++ b/estimagic/optimization/optimize.py
@@ -7,10 +7,12 @@ from multiprocessing import Queue
 from pathlib import Path
 
 import numpy as np
+import pandas as pd
 import pygmo as pg
 from scipy.optimize import minimize as scipy_minimize
 
 from estimagic.dashboard.server_functions import run_server
+from estimagic.differentiation.differentiation import gradient
 from estimagic.optimization.process_constraints import process_constraints
 from estimagic.optimization.reparametrize import reparametrize_from_internal
 from estimagic.optimization.reparametrize import reparametrize_to_internal
@@ -147,7 +149,17 @@ def minimize(
     params = _process_params(params)
     fitness_eval = criterion(params, *criterion_args, **criterion_kwargs)
     constraints = process_constraints(constraints, params)
-    internal_params = reparametrize_to_internal(params, constraints)
+    internal_params = reparametrize_to_internal(params, constraints, None)
+
+    scaling_factor = calculate_scaling_factor(
+        criterion,
+        params,
+        internal_params,
+        constraints,
+        general_options,
+        criterion_args,
+        criterion_kwargs,
+    )
 
     queue = Queue() if dashboard else None
     if dashboard:
@@ -172,6 +184,7 @@ def minimize(
         params=params,
         internal_params=internal_params,
         constraints=constraints,
+        scaling_factor=scaling_factor,
         algorithm=algorithm,
         algo_options=algo_options,
         general_options=general_options,
@@ -191,6 +204,7 @@ def _minimize(
     params,
     internal_params,
     constraints,
+    scaling_factor,
     algorithm,
     algo_options,
     general_options,
@@ -233,11 +247,12 @@ def _minimize(
             the updated parameter Series will be supplied later.
 
     """
-    internal_criterion = _create_internal_criterion(
+    internal_criterion = create_internal_criterion(
         criterion=criterion,
         params=params,
         internal_params=internal_params,
         constraints=constraints,
+        scaling_factor=scaling_factor,
         criterion_args=criterion_args,
         criterion_kwargs=criterion_kwargs,
         queue=queue,
@@ -263,10 +278,12 @@ def _minimize(
         algo = _create_algorithm(algo_name, algo_options, origin)
         pop = _create_population(prob, algo_options, internal_params)
         evolved = algo.evolve(pop)
-        result = _process_results(evolved, params, internal_params, constraints, origin)
+        result = _process_results(
+            evolved, params, internal_params, constraints, origin, scaling_factor
+        )
     elif origin == "scipy":
         bounds = _get_scipy_bounds(internal_params)
-        x0 = _x_from_params(params, constraints)
+        x0 = _x_from_params(params, constraints, scaling_factor)
         minimized = scipy_minimize(
             internal_criterion,
             x0,
@@ -275,7 +292,7 @@ def _minimize(
             options=algo_options,
         )
         result = _process_results(
-            minimized, params, internal_params, constraints, origin
+            minimized, params, internal_params, constraints, origin, scaling_factor
         )
     else:
         raise ValueError("Invalid algorithm requested.")
@@ -283,11 +300,12 @@ def _minimize(
     return result
 
 
-def _create_internal_criterion(
+def create_internal_criterion(
     criterion,
     params,
     internal_params,
     constraints,
+    scaling_factor,
     criterion_args,
     criterion_kwargs,
     queue,
@@ -295,7 +313,7 @@ def _create_internal_criterion(
     c = np.ones(1, dtype=int)
 
     def internal_criterion(x, counter=c):
-        p = _params_from_x(x, internal_params, constraints, params)
+        p = _params_from_x(x, internal_params, constraints, params, scaling_factor)
         fitness_eval = criterion(p, *criterion_args, **criterion_kwargs)
         if queue is not None:
             queue.put(QueueEntry(iteration=counter[0], params=p, fitness=fitness_eval))
@@ -305,15 +323,23 @@ def _create_internal_criterion(
     return internal_criterion
 
 
-def _params_from_x(x, internal_params, constraints, params):
+def _params_from_x(x, internal_params, constraints, params, scaling_factor):
     internal_params = internal_params.copy(deep=True)
-    internal_params["value"] = x
-    updated_params = reparametrize_from_internal(internal_params, constraints, params)
+    # :func:`internal_criterion` always assumes that `x` is a NumPy array, but if we
+    # pass the internal criterion function to :func:`gradient`, x is a DataFrame.
+    # Setting a series to a DataFrame will convert the column "value" to object type
+    # which causes trouble in following NumPy routines assuming a numeric type.
+    internal_params["value"] = x["value"] if isinstance(x, pd.DataFrame) else x
+    updated_params = reparametrize_from_internal(
+        internal_params, constraints, params, scaling_factor
+    )
     return updated_params
 
 
-def _x_from_params(params, constraints):
-    return reparametrize_to_internal(params, constraints)["value"].to_numpy()
+def _x_from_params(params, constraints, scaling_factor):
+    return reparametrize_to_internal(params, constraints, scaling_factor)[
+        "value"
+    ].to_numpy()
 
 
 def _process_params(params):
@@ -395,7 +421,7 @@ def _create_population(problem, algo_options, internal_params):
     return pop
 
 
-def _process_results(res, params, internal_params, constraints, origin):
+def _process_results(res, params, internal_params, constraints, origin, scaling_factor):
     """Convert optimization results into json serializable dictionary.
     Args:
         res: Result from numerical optimizer.
@@ -414,7 +440,85 @@ def _process_results(res, params, internal_params, constraints, origin):
     elif origin in ["pygmo", "nlopt"]:
         x = res.champion_x
         res_dict = {"fun": res.champion_f[0]}
-    params = _params_from_x(x, internal_params, constraints, params)
+    params = _params_from_x(x, internal_params, constraints, params, scaling_factor)
 
     res_dict["internal_x"] = x.tolist()
     return res_dict, params
+
+
+def calculate_scaling_factor(
+    criterion,
+    params,
+    internal,
+    constraints,
+    general_options,
+    criterion_args,
+    criterion_kwargs,
+):
+    """Calculate the scaling factor for the internal parameters.
+
+    There are multiple ways the user is able to rescale the parameter vector which is
+    specified in ``general_options["scaling"]``.
+
+    - ``None`` (default): No scaling happens.
+    - ``"start_values"``: Divide parameters which are not in ``[-1, 1]`` by their
+      starting values.
+    - ``"gradient"``: Divide parameters which are not in ``[-1e-2, 1e-2]`` by the
+      inverse of the gradient. By default, the computation method is ``"central"`` with
+      extrapolation set to ``False``. The user can change the defaults by passing
+      ``scaling_gradient_method`` or ``scaling_gradient_extrapolation``. See
+      :func:`~estimagic.differentiation.differentiation.gradient` for more details.
+
+    Note that the scaling factor should be defined such that unscaling is done by
+    multiplying the scaling factor. This simplifies :func:`_rescale_from_internal`.
+
+    Args:
+        criterion (func): Criterion function.
+        params (DataFrame): See :ref:`params`.
+        internal (DataFrame): See :ref:`params`.
+        constraints (dict): Dictionary containing constraints.
+        general_options (dict): General options. See :ref:`estimation_general_options`.
+        criterion_args (list): List of arguments of the criterion function.
+        crtierion_kwargs (dict): Dictionary of keyword arguments of the criterion
+            function.
+
+    Returns:
+        scaling_factor (Series): Scaling Factor.
+
+    """
+    scaling = general_options.get("scaling", None)
+
+    if scaling is None:
+        scaling_factor = np.ones(len(internal))
+    elif scaling == "start_values":
+        scaling_factor = internal["value"].abs().clip(1)
+    elif scaling == "gradient":
+        method = general_options.get("scaling_gradient_method", "central")
+        extrapolation = general_options.get("scaling_gradient_extrapolation", False)
+
+        internal_criterion = create_internal_criterion(
+            criterion,
+            params,
+            internal,
+            constraints,
+            None,
+            criterion_args,
+            criterion_kwargs,
+            queue=None,
+        )
+
+        gradients = gradient(
+            internal_criterion,
+            internal,
+            method,
+            extrapolation,
+            criterion_args,
+            criterion_kwargs,
+        )
+        scaling_factor = np.clip(1 / gradients.abs(), 1e-2, None)
+    else:
+        raise NotImplementedError(f"Scaling method {scaling} is not implemented.")
+
+    scaling_factor = pd.Series(scaling_factor, index=internal.index)
+
+    return scaling_factor
diff --git a/estimagic/optimization/reparametrize.py b/estimagic/optimization/reparametrize.py
index f005d78..9a33111 100644
--- a/estimagic/optimization/reparametrize.py
+++ b/estimagic/optimization/reparametrize.py
@@ -12,7 +12,7 @@ from estimagic.optimization.utilities import robust_cholesky
 from estimagic.optimization.utilities import sdcorr_params_to_matrix
 
 
-def reparametrize_to_internal(params, constraints):
+def reparametrize_to_internal(params, constraints, scaling_factor):
     """Convert a params DataFrame to an internal_params DataFrame.
 
     The internal params df is shorter because it does not contain fixed parameters.
@@ -61,6 +61,8 @@ def reparametrize_to_internal(params, constraints):
     internal = internal.loc[~(internal["_fixed"])].copy(deep=True)
     internal.drop(columns="_fixed", axis=1, inplace=True)
 
+    internal = _scaling_to_internal(internal, scaling_factor)
+
     invalid = internal.query("lower >= upper | lower > value | upper < value")
     assert (
         len(invalid) == 0
@@ -69,7 +71,9 @@ def reparametrize_to_internal(params, constraints):
     return internal
 
 
-def reparametrize_from_internal(internal_params, constraints, original_params):
+def reparametrize_from_internal(
+    internal_params, constraints, original_params, scaling_factor
+):
     """Convert an internal_params DataFrame to a Series with valid parameters.
 
     The parameter values are constructed from the 'value' column of internal_params.
@@ -87,7 +91,11 @@ def reparametrize_from_internal(internal_params, constraints, original_params):
 
     """
     external = original_params.copy(deep=True)
+
+    internal_params = _scaling_from_internal(internal_params, scaling_factor)
+
     external.update(internal_params["value"])
+
     external["_fixed"] = True
     external.loc[internal_params.index, "_fixed"] = False
 
@@ -420,3 +428,23 @@ def _fixed_from_internal(params_subset, constr):
     if value is not None:
         res["value"] = value
     return res["value"]
+
+
+def _scaling_to_internal(internal, scaling_factor):
+    """Scale parameters by division with the scaling factor."""
+    if scaling_factor is not None:
+        internal[["value", "lower", "upper"]] = internal[
+            ["value", "lower", "upper"]
+        ].divide(scaling_factor, axis="index")
+
+    return internal
+
+
+def _scaling_from_internal(internal, scaling_factor):
+    """Unscale parameters by multiplication with the scaling factor."""
+    if scaling_factor is not None:
+        internal[["value", "lower", "upper"]] = internal[
+            ["value", "lower", "upper"]
+        ].multiply(scaling_factor, axis="index")
+
+    return internal
diff --git a/setup.py b/setup.py
index 06d9410..52a2736 100644
--- a/setup.py
+++ b/setup.py
@@ -1,7 +1,7 @@
 from setuptools import find_packages
 from setuptools import setup
 
-version = "0.0.14"
+version = "0.0.15"
 
 setup(
     name="estimagic",
diff --git a/tox.ini b/tox.ini
index 69beb8d..d8a92bc 100644
--- a/tox.ini
+++ b/tox.ini
@@ -28,7 +28,6 @@ commands = pytest estimagic
 passenv = USERPROFILE SSH_AUTH_SOCK
 deps =
     pre-commit
-    doc8
     restructuredtext_lint
     Pygments
 commands =
@@ -70,3 +69,4 @@ norecursedirs =
 filterwarnings =
     ignore: Using or importing the ABCs from 'collections'
     ignore: the imp module is deprecated
+    ignore: indexing past lexsort depth may impact performance.


Test patch:
diff --git a/estimagic/tests/optimization/test_reparametrize.py b/estimagic/tests/optimization/test_reparametrize.py
index 35cb103..e4515d0 100644
--- a/estimagic/tests/optimization/test_reparametrize.py
+++ b/estimagic/tests/optimization/test_reparametrize.py
@@ -1,4 +1,3 @@
-import warnings
 from pathlib import Path
 
 import numpy as np
@@ -76,14 +75,11 @@ for ext, int_ in zip(external, internal):
 def test_reparametrize_to_internal(params, expected_internal, category):
     constr = constraints(params)
     cols = ["value", "lower", "upper"]
-    with warnings.catch_warnings():
-        warnings.filterwarnings(
-            "ignore", message="indexing past lexsort depth may impact performance."
-        )
-        calculated = reparametrize_to_internal(params, constr)
-        assert_frame_equal(
-            calculated.loc[category, cols], expected_internal.loc[category, cols]
-        )
+
+    calculated = reparametrize_to_internal(params, constr, None)
+    assert_frame_equal(
+        calculated.loc[category, cols], expected_internal.loc[category, cols]
+    )
 
 
 to_test = []
@@ -95,13 +91,8 @@ for int_, ext in zip(internal, external):
 @pytest.mark.parametrize("internal, expected_external, category", to_test)
 def test_reparametrize_from_internal(internal, expected_external, category):
     constr = constraints(expected_external)
-    with warnings.catch_warnings():
-        warnings.filterwarnings(
-            "ignore", message="indexing past lexsort depth may impact performance."
-        )
-        calculated = reparametrize_from_internal(internal, constr, expected_external)[
-            "value"
-        ]
-        assert_series_equal(
-            calculated[category], expected_external.loc[category, "value"]
-        )
+
+    calculated = reparametrize_from_internal(internal, constr, expected_external, None)[
+        "value"
+    ]
+    assert_series_equal(calculated[category], expected_external.loc[category, "value"])