Fix issues related to providing a Jacobian function (#973)

* Updated the Minimizer.least_squares method to wrap the Jacobian function
- The Jacobian function provided by the user is now handled similarly to how it is handled in, e.g., Minimizer.leastsq.
- If the Jacobian function is called 'Dfun' in the keyword arguments, then it is also accepted (Minimizer.scalar_minimize already does this).

* Updated the Minimizer.__jacobian method signature to fix exception
- The signature now matches the signature of Minimizer.__residual. This was changed because the following exception was being raised whenever the user provided a function for calculating the Jacobian: "TypeError: Minimizer.__jacobian() got an unexpected keyword argument 'apply_bounds_transformation'".

* Updated the Minimizer.__jacobian method to use the new parameter
- The updating of the parameter values is now similar to how it is done in Minimizer.__residual.

* Updated test_least_squares_jacobian_types to handle the wrapped Jacobian function
- The params argument of the jac_array function is now an actual Parameters object like it is in the f function, which means that the parameters can be accessed via string keys.

* Updated the coerce_float64 function
- Added explicit support for sparse matrices (e.g., Block Sparse Row) and LinearOperator objects. This change is needed because the updated test_least_squares_jacobian_types would otherwise fail. Similar to how these types are handled when calculating the Hessian matrix in Minimizer.least_squares.

* Renamed the Minimizer.__jacobian method to avoid a bug related to pickling
- Attempting to make use of multiprocessing.Pool while providing, e.g., Minimizer.least_squares a function for calculating the Jacobian did not work prior to the renaming of the method due to an issue related to pickling and unpickling of private methods. See the (currently) open issue (python/cpython#77188) and pull request (python/cpython#21480) for more information.

* Added a test for when the "least_squares" method is provided a function for the Jacobian
- Verifies that a) the fitted parameters have the expected values, b) the jac function is actually called, and c) that there are indeed fewer function evaluations as a result.

* Updated the docstring of Minimizer._jacobian
- The docstring is updated to include the parameters and return values.

* Added tests for (un)pickling Minimizer when a Jacobian function is provided

* Updated timeout value in test_jacobian_with_forkingpickler
- Increased the timeout value so that slow runners can actually succeed instead of timing out. If no timeout was specified, then a failing test could end up hanging instead.

* Added an example that demonstrates the benefits of providing a Jacobian function

* Updated the list of contributors

AUTHORS.txt

@@ -55,7 +55,7 @@ Roam, Alexander Stark, Alexandre Beelen, Andrey Aristov, Nicholas Zobrist,
 Ethan Welty, Julius Zimmermann, Mark Dean, Arun Persaud, Ray Osborn, @lneuhaus,
 Marcel Stimberg, Yoshiera Huang, Leon Foks, Sebastian Weigand, Florian LB,
 Michael Hudson-Doyle, Ruben Verweij, @jedzill4, @spalato, Jens Hedegaard Nielsen,
-Martin Majli, Kristian Meyer, @azelcer, Ivan Usov, and many others.
+Martin Majli, Kristian Meyer, @azelcer, Ivan Usov, Ville Yrjänä, and many others.
 
 The lmfit code obviously depends on, and owes a very large debt to the code
 in scipy.optimize. Several discussions on the SciPy-user and lmfit mailing

examples/example_fit_jacobian_benchmark.py

@@ -0,0 +1,178 @@
+"""
+Benchmarks of methods with and without computing the Jacobian analytically
+==========================================================================
+
+Providing a function that calculates the Jacobian matrix analytically can
+reduce the time spent finding a solution. The results from benchmarks comparing
+two methods (``leastsq`` and ``least_squares``) with and without a function to
+calculate the Jacobian matrix analytically are presented below.
+
+First we define the model function, the residual function, and the appropriate
+Jacobian functions:
+"""
+from timeit import timeit
+from types import SimpleNamespace
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from lmfit import Parameters, minimize
+
+NUM_JACOBIAN_CALLS = 0
+
+
+def func(var, x):
+    return var[0] * np.exp(-var[1]*x) + var[2]
+
+
+def residual(pars, x, data):
+    a, b, c = pars['a'], pars['b'], pars['c']
+    model = func((a, b, c), x)
+    return model - data
+
+
+def dfunc(pars, x, data):
+    global NUM_JACOBIAN_CALLS
+    NUM_JACOBIAN_CALLS += 1
+
+    a, b = pars['a'], pars['b']
+    v = np.exp(-b*x)
+    return np.array([v, -a*x*v, np.ones(len(x))])
+
+
+def jacfunc(pars, x, data):
+    global NUM_JACOBIAN_CALLS
+    NUM_JACOBIAN_CALLS += 1
+
+    a, b = pars['a'], pars['b']
+    v = np.exp(-b*x)
+    jac = np.ones((len(x), 3), dtype=np.float64)
+    jac[:, 0] = v
+    jac[:, 1] = -a * x * v
+    return jac
+
+
+a, b, c = 2.5, 1.3, 0.8
+
+x = np.linspace(0, 4, 50)
+y = func([a, b, c], x)
+
+data = y + 0.15*np.random.RandomState(seed=2021).normal(size=x.size)
+
+
+###############################################################################
+# Then we define the different cases to benchmark (i.e., different methods with
+# and without a function to calculate the Jacobian analytically) and the number
+# of repetitions per case:
+cases = (
+    dict(
+        method='leastsq',
+    ),
+    dict(
+        method='leastsq',
+        Dfun=dfunc,
+        col_deriv=1,
+    ),
+    dict(
+        method='least_squares',
+    ),
+    dict(
+        method='least_squares',
+        jac=jacfunc,
+    ),
+)
+
+num_repeats = 100
+results = []
+
+for kwargs in cases:
+    params = Parameters()
+    params.add('a', value=10)
+    params.add('b', value=10)
+    params.add('c', value=10)
+
+    wrapper = lambda: minimize(
+        residual,
+        params,
+        args=(x,),
+        kws={'data': data},
+        **kwargs,
+    )
+    time = timeit(wrapper, number=num_repeats) / num_repeats
+
+    NUM_JACOBIAN_CALLS = 0
+    fit = wrapper()
+
+    results.append(SimpleNamespace(
+        time=time,
+        num_jacobian_calls=NUM_JACOBIAN_CALLS,
+        fit=fit,
+        kwargs=kwargs,
+    ))
+
+
+###############################################################################
+# Finally, we present the results:
+labels = []
+
+for result in results:
+    label = result.kwargs['method']
+    if result.num_jacobian_calls > 0:
+        label += ' with Jac.'
+
+    labels.append(label)
+
+label_width = max(map(len, labels))
+lines = [
+    '| '
+    + ' | '.join([
+        'Method'.ljust(label_width),
+        'Avg. time (ms)',
+        '# func. (+ Jac.) calls',
+        'Chi-squared',
+        'a'.ljust(5),
+        'b'.ljust(5),
+        'c'.ljust(6),
+    ])
+    + '|'
+]
+
+print(f'The "true" parameters are: a = {a:.3f}, b = {b:.3f}, c = {c:.3f}\n')
+fig, ax = plt.subplots()
+ax.plot(x, data, marker='.', linestyle='none', label='data')
+
+for (result, label) in zip(results, labels):
+    linestyle = '-'
+    if result.num_jacobian_calls > 0:
+        linestyle = '--'
+
+    a = result.fit.params['a'].value
+    b = result.fit.params['b'].value
+    c = result.fit.params['c'].value
+    y = func([a, b, c], x)
+    ax.plot(x, y, label=label, alpha=0.5, linestyle=linestyle)
+
+    columns = [
+        label.ljust(label_width),
+        f'{result.time * 1000:.2f}'.ljust(14),
+        (
+            f'{result.fit.nfev}'
+            + (
+                f' (+{result.num_jacobian_calls})'
+                if result.num_jacobian_calls > 0 else
+                ''
+            )
+        ).ljust(22),
+        f'{result.fit.chisqr:.3f}'.ljust(11),
+        f'{a:.3f}'.ljust(5, '0'),
+        f'{b:.3f}'.ljust(5, '0'),
+        f'{c:.3f}'.ljust(5, '0'),
+    ]
+    lines.append('| ' + ' | '.join(columns) + ' |')
+
+lines.insert(1, '|-' + '-|-'.join('-' * len(col) for col in columns) + '-|')
+print('\n'.join(lines))
+
+ax.set_xlabel('x')
+ax.set_ylabel('y')
+ax.legend()

lmfit/minimizer.py

@@ -553,20 +553,37 @@ def __residual(self, fvars, apply_bounds_transformation=True):
         else:
             return coerce_float64(out, nan_policy=self.nan_policy)
 
-    def __jacobian(self, fvars):
+    def _jacobian(self, fvars, apply_bounds_transformation=True):
         """Return analytical jacobian to be used with Levenberg-Marquardt.
 
         modified 02-01-2012 by Glenn Jones, Aberystwyth University
         modified 06-29-2015 by M Newville to apply gradient scaling for
         bounded variables (thanks to JJ Helmus, N Mayorov)
 
+        Parameters
+        ----------
+        fvars : numpy.ndarray
+            Array of new parameter values suggested by the minimizer.
+        apply_bounds_transformation : bool, optional
+            Whether to apply lmfits parameter transformation to constrain
+            parameters (default is True). This is needed for solvers
+            without built-in support for bounds.
+
+        Returns
+        -------
+        numpy.ndarray
+             The evaluated Jacobian matrix for given `fvars`.
+
         """
         pars = self.result.params
         grad_scale = np.ones_like(fvars)
         for ivar, name in enumerate(self.result.var_names):
             val = fvars[ivar]
-            pars[name].value = pars[name].from_internal(val)
-            grad_scale[ivar] = pars[name].scale_gradient(val)
+            if apply_bounds_transformation:
+                pars[name].value = pars[name].from_internal(val)
+                grad_scale[ivar] = pars[name].scale_gradient(val)
+            else:
+                pars[name].value = val
 
         pars.update_constraints()
 
@@ -931,7 +948,7 @@ def scalar_minimize(self, method='Nelder-Mead', params=None, max_nfev=None,
         # Wrap Jacobian function to deal with bounds
         if 'jac' in fmin_kws:
             self.jacfcn = fmin_kws.pop('jac')
-            fmin_kws['jac'] = self.__jacobian
+            fmin_kws['jac'] = self._jacobian
 
         # Ignore jac argument for methods that do not support it
         if 'jac' in fmin_kws and method not in ('CG', 'BFGS', 'Newton-CG',
@@ -1532,6 +1549,13 @@ def least_squares(self, params=None, max_nfev=None, **kws):
         least_squares_kws.update(self.kws)
         least_squares_kws.update(kws)
 
+        if least_squares_kws.get('Dfun', None) is not None:
+            least_squares_kws['jac'] = least_squares_kws.pop('Dfun')
+
+        if callable(least_squares_kws['jac']):
+            self.jacfcn = least_squares_kws['jac']
+            least_squares_kws['jac'] = self._jacobian
+
         least_squares_kws['kwargs'].update({'apply_bounds_transformation': False})
         result.call_kws = least_squares_kws
 
@@ -1639,7 +1663,7 @@ def leastsq(self, params=None, max_nfev=None, **kws):
         if lskws['Dfun'] is not None:
             self.jacfcn = lskws['Dfun']
             self.col_deriv = lskws['col_deriv']
-            lskws['Dfun'] = self.__jacobian
+            lskws['Dfun'] = self._jacobian
 
         # suppress runtime warnings during fit and error analysis
         orig_warn_settings = np.geterr()
@@ -2386,7 +2410,12 @@ def coerce_float64(arr, nan_policy='raise', handle_inf=True,
     lists of numbers, pandas.Series, h5py.Datasets, and many other array-like
     Python objects
     """
-    if np.iscomplexobj(arr):
+    if issparse(arr):
+        arr = arr.toarray().astype(np.float64)
+    elif isinstance(arr, LinearOperator):
+        identity = np.eye(arr.shape[1], dtype=np.float64)
+        arr = (arr * identity).astype(np.float64)
+    elif np.iscomplexobj(arr):
         arr = np.asarray(arr, dtype=np.complex128).view(np.float64)
     else:
         arr = np.asarray(arr, dtype=np.float64)