diff --git a/AUTHORS.txt b/AUTHORS.txt
index 7ad1d47e..889c9331 100644
--- a/AUTHORS.txt
+++ b/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
diff --git a/examples/example_fit_jacobian_benchmark.py b/examples/example_fit_jacobian_benchmark.py
new file mode 100644
index 00000000..9fdcca0b
--- /dev/null
+++ b/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()
diff --git a/lmfit/minimizer.py b/lmfit/minimizer.py
index eb3f8e3b..fd20afe2 100644
--- a/lmfit/minimizer.py
+++ b/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)
diff --git a/tests/test_jacobian_pickling.py b/tests/test_jacobian_pickling.py
new file mode 100644
index 00000000..ff5d9e33
--- /dev/null
+++ b/tests/test_jacobian_pickling.py
@@ -0,0 +1,129 @@
+"""Tests for (un)pickling the Minimizer class when providing a Jacobian function."""
+from multiprocessing import get_context
+from pickle import dumps, loads
+
+import numpy as np
+
+from lmfit import Minimizer, Parameters, minimize
+
+
+def func(pars, x, data=None):
+    a, b, c = pars['a'], pars['b'], pars['c']
+    model = a * np.exp(-b*x) + c
+    if data is None:
+        return model
+    return model - data
+
+
+def dfunc(pars, x, data=None):
+    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=None):
+    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
+
+
+def f(var, x):
+    return var[0] * np.exp(-var[1]*x) + var[2]
+
+
+def wrapper(args):
+    (
+        method,
+        x,
+        data,
+    ) = args
+    params = Parameters()
+    params.add('a', value=10)
+    params.add('b', value=10)
+    params.add('c', value=10)
+
+    return minimize(
+        func,
+        params,
+        method=method,
+        args=(x,),
+        kws={'data': data},
+        **(
+            {'Dfun': dfunc, 'col_deriv': 1}
+            if method == 'leastsq' else
+            {'jac': jacfunc}
+        ),
+    )
+
+
+def test_jacobian_with_pickle():
+    """Test using pickle.dumps/loads."""
+    params = Parameters()
+    params.add('a', value=10)
+    params.add('b', value=10)
+    params.add('c', value=10)
+
+    a, b, c = 2.5, 1.3, 0.8
+    x = np.linspace(0, 4, 50)
+    y = f([a, b, c], x)
+    np.random.seed(2021)
+    data = y + 0.15*np.random.normal(size=x.size)
+
+    for method in ('leastsq', 'least_squares'):
+        if method == 'leastsq':
+            kwargs = {'Dfun': dfunc, 'col_deriv': 1}
+        else:
+            kwargs = {'jac': jacfunc}
+
+        min = Minimizer(
+            func,
+            params,
+            fcn_args=(x,),
+            fcn_kws={'data': data},
+            **kwargs,
+        )
+
+        pickled = dumps(min)
+        unpickled = loads(pickled)
+
+        out = unpickled.minimize(method=method)
+
+        assert np.isclose(out.params['a'], 2.5635, atol=1e-4)
+        assert np.isclose(out.params['b'], 1.3585, atol=1e-4)
+        assert np.isclose(out.params['c'], 0.8241, atol=1e-4)
+
+
+def test_jacobian_with_forkingpickler():
+    """Test using multiprocessing.Pool, which uses a subclass of pickle.Pickler."""
+    a, b, c = 2.5, 1.3, 0.8
+    x = np.linspace(0, 4, 50)
+    y = f([a, b, c], x)
+    np.random.seed(2021)
+    data = y + 0.15*np.random.normal(size=x.size)
+
+    with get_context(method='spawn').Pool(1) as pool:
+        iterator = pool.imap_unordered(
+            wrapper,
+            (
+                (
+                    method,
+                    x,
+                    data,
+                )
+                for method in ('leastsq', 'least_squares')
+            ),
+            chunksize=1,
+        )
+
+        while True:
+            try:
+                out = iterator.next(timeout=30)
+            except StopIteration:
+                break
+
+            assert np.isclose(out.params['a'], 2.5635, atol=1e-4)
+            assert np.isclose(out.params['b'], 1.3585, atol=1e-4)
+            assert np.isclose(out.params['c'], 0.8241, atol=1e-4)
diff --git a/tests/test_least_squares.py b/tests/test_least_squares.py
index bb592c72..301a75c0 100644
--- a/tests/test_least_squares.py
+++ b/tests/test_least_squares.py
@@ -138,7 +138,7 @@ def f(params):
     # numpy.ndarray, scipy.sparse.spmatrix, scipy.sparse.linalg.LinearOperator
     # J = [ 2x - 2a , 2y - 2b ]
     def jac_array(params, *args, **kwargs):
-        return np.column_stack((2 * params[0] - 2 * a, 2 * params[1] - 2 * b))
+        return np.column_stack((2 * params['x'] - 2 * a, 2 * params['y'] - 2 * b))
 
     def jac_sparse(params, *args, **kwargs):
         return bsr_matrix(jac_array(params, *args, **kwargs))
@@ -169,3 +169,47 @@ def jac_operator(params, *args, **kwargs):
     assert_allclose(result.covar, result_array.covar)
     assert_allclose(result.covar, result_sparse.covar)
     assert_allclose(result.covar, result_operator.covar)
+
+
+def test_least_squares_jacobian():
+    pars = lmfit.Parameters()
+    pars.add('x0', value=2.0)
+    pars.add('x1', value=2.0, min=1.5)
+
+    global jac_count
+
+    jac_count = 0
+
+    def resid(params):
+        x0 = params['x0']
+        x1 = params['x1']
+        return np.array([10 * (x1 - x0*x0), 1-x0])
+
+    def jac(params):
+        global jac_count
+        jac_count += 1
+        x0 = params['x0']
+        return np.array([[-20*x0, 10], [-1, 0]])
+
+    out0 = lmfit.minimize(resid, pars, method='least_squares', jac='2-point')
+
+    assert_allclose(out0.params['x0'], 1.2243, rtol=1.0e-4)
+    assert_allclose(out0.params['x1'], 1.5000, rtol=1.0e-4)
+    assert jac_count == 0
+
+    out1 = lmfit.minimize(resid, pars, method='least_squares', jac=jac)
+
+    assert_allclose(out1.params['x0'], 1.2243, rtol=1.0e-4)
+    assert_allclose(out1.params['x1'], 1.5000, rtol=1.0e-4)
+    assert out1.nfev < out0.nfev
+    assert jac_count > 5
+
+    jac_count = 0
+
+    out2 = lmfit.minimize(resid, pars, method='least_squares', Dfun=jac)
+
+    assert_allclose(out2.params['x0'], 1.2243, rtol=1.0e-4)
+    assert_allclose(out2.params['x1'], 1.5000, rtol=1.0e-4)
+    assert out2.nfev < out0.nfev
+    assert out2.nfev == out1.nfev
+    assert jac_count > 5