codecov[bot]: # [Codecov](https://codecov.io/gh/lmfit/lmfit-py/pull/589?src=pr&el=h1) Report
> Merging [#589](https://codecov.io/gh/lmfit/lmfit-py/pull/589?src=pr&el=desc) into [master](https://codecov.io/gh/lmfit/lmfit-py/commit/d7b408baa64821f4356cff9c32a199bc53becd93?src=pr&el=desc) will **increase** coverage by `0.03%`.
> The diff coverage is `100%`.

[![Impacted file tree graph](https://codecov.io/gh/lmfit/lmfit-py/pull/589/graphs/tree.svg?width=650&token=hvITDAARAt&height=150&src=pr)](https://codecov.io/gh/lmfit/lmfit-py/pull/589?src=pr&el=tree)

```diff
@@            Coverage Diff            @@
##           master    #589      +/-   ##
=========================================
+ Coverage   90.46%   90.5%   +0.03%     
=========================================
  Files          11      11              
  Lines        3263    3274      +11     
=========================================
+ Hits         2952    2963      +11     
  Misses        311     311
```


| [Impacted Files](https://codecov.io/gh/lmfit/lmfit-py/pull/589?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [lmfit/minimizer.py](https://codecov.io/gh/lmfit/lmfit-py/pull/589/diff?src=pr&el=tree#diff-bG1maXQvbWluaW1pemVyLnB5) | `91.59% <100%> (+0.11%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/lmfit/lmfit-py/pull/589?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/lmfit/lmfit-py/pull/589?src=pr&el=footer). Last update [d7b408b...36a022f](https://codecov.io/gh/lmfit/lmfit-py/pull/589?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

newville: @ezwelty  Like @reneeotten I'm mostly OK with this. 

I have two concerns:  First, I think that @reneeotten are unlikely to be able to support this.   Adding unsupported features seems like maybe not such a great idea.   Are you willing to support this?

The example you gave is really close to trivial.  That's OK for a test -- it shows the code "works".    But it leads to a second concern:  What happens if someone places bounds on a parameter or even fixes a parameter?    Like, why is the sparse jacobian a 2-d array?  Why is not more naturally (for lmfit) a dictionary with parameter names as keys and 1d arrays of the sparsity of that derivative component?

reneeotten: @newville I am not sure I understand your concerns completely. To me it seems the "only" thing this PR does is making sure that whatever type of Jacobian one uses, it will return the Hessian matrix in the same way. 

The ```least_squares``` method already supports bounds and so the resulting Hessian will be already correct (right now we also don't transform the ```result.jac``` from ```least_squares```). In the test-case it is a 2x2 matrix because there are two fitting parameters, if you'd fix one of them it should be only one value (but I agree it wouldn't hurt to add another test where parameters have bounds and some are fixed, just to make sure).

The resulting Hessian matrix will be used in the "lmfit-way" to calculate the uncertainties using:
```
result.covar = np.linalg.inv(hess)
self._calculate_uncertainties_correlations()
```
and I don't see anything wrong with that; but perhaps, I'm missing something here and, if so, can you elaborate a bit?
newville: @reneeotten   I agree that bounds ought to be OK.  Maybe that would be good to test. 

If you have a Parameter with `vary=False` you now have fewer variables than Parameters.  The order of the variables is predictable, but the Jacobian (and so its sparsity) has changed size.    If the Parameter is fixed, the Jacobin is not just sparse, it is non-existent.

With this PR (and #588), the user can provide a `jac_sparsity` as input.   That is OK.  But that is a 2D array that is Nvariables X Nobservations.

Nothing else in `lmfit` accesses variables by index, and it can be confusing which variable is in "position 2" -- the user can see this after the fit from `var_names`.   But they don't have to know this order when they do the fit. 

If the user instead provided a dictionary of 1D derivatives and sparsities or made these as Parameter attributes then `minimize()` could construct the internal `jac_sparsity` 2D array that gets sent to `trf`.  

Does that make sense?  






ezwelty: @newville I see your point.

With this PR, the user can now pass `jac_sparsity` and all allowed flavors of `jac` to `scipy.optimize.least_squares` via optional keyword arguments. It is, however, also up to the user to make sure that these are consistent with the order and inclusion (`vary=True/False`) of parameters passed by `lmfit`. (Bounds do not change the dimensions of the Jacobian). This isn't only true of this PR – `jac` (returning an array) is implicitly supported previous to this PR, as well as passing `x_scale` or `diff_step`, and the same disclaimer holds. I would suggest that protecting the user from passing misshapen optional arguments to `least_squares` is beyond the scope of this PR, and this PR will be needed regardless. In the meantime, I think it is more helpful to empower users to make use of these optional features than to block them from doing so.

If folks agree, I can brush up the pull request as requested by @reneeotten.
newville: @ezwelty @reneeotten Any thoughts on the status of this PR?  I'm OK with merging, if the requested changes and conflicts are resolved.  

It does sort of seem like we should improve the interface for Jacobians, but I think that is my job and shouldn't hold up this PR. 


ezwelty: @newville Sorry, I wasn't sure I was cleared to move forward. I will brush up the PR as requested by @reneeotten.
diff --git a/lmfit/minimizer.py b/lmfit/minimizer.py
index ce2cbed5..22c275be 100644
--- a/lmfit/minimizer.py
+++ b/lmfit/minimizer.py
@@ -29,6 +29,8 @@ from scipy.optimize import dual_annealing as scipy_dual_annealing
 from scipy.optimize import leastsq as scipy_leastsq
 from scipy.optimize import minimize as scipy_minimize
 from scipy.optimize import shgo as scipy_shgo
+from scipy.sparse import issparse
+from scipy.sparse.linalg import LinearOperator
 from scipy.stats import cauchy as cauchy_dist
 from scipy.stats import norm as norm_dist
 from scipy.version import version as scipy_version
@@ -1408,7 +1410,21 @@ class Minimizer(object):
 
             # calculate the cov_x and estimate uncertainties/correlations
             try:
-                hess = np.matmul(ret.jac.T, ret.jac)
+                if issparse(ret.jac):
+                    hess = (ret.jac.T * ret.jac).toarray()
+                elif isinstance(ret.jac, LinearOperator):
+                    identity = np.eye(ret.jac.shape[1], dtype=ret.jac.dtype)
+                    # TODO: Remove try-except when scipy < 1.4.0 support dropped
+                    try:
+                        # For scipy >= 1.4.0 (with Linear Operator transpose)
+                        # https://github.com/scipy/scipy/pull/9064
+                        hess = (ret.jac.T * ret.jac) * identity
+                    except AttributeError:
+                        # For scipy < 1.4.0 (without Linear Operator transpose)
+                        jac = ret.jac * identity
+                        hess = np.matmul(jac.T, jac)
+                else:
+                    hess = np.matmul(ret.jac.T, ret.jac)
                 result.covar = np.linalg.inv(hess)
                 self._calculate_uncertainties_correlations()
             except LinAlgError:

Passing jac_sparsity to least_squares breaks minimize()
#### Description

Calls of the form:

```python
lmfit.minimize(..., method="least_squares", jac_sparsity=obj)
```

result in the Jacobian being returned as a sparse matrix, which breaks the calculation of the Hessian:

https://github.com/lmfit/lmfit-py/blob/d7b408baa64821f4356cff9c32a199bc53becd93/lmfit/minimizer.py#L1430

This should be replaced with something like:

```python
if scipy.sparse.issparse(ret.jac):
    jac = ret.jac.todense()
else:
    jac = ret.jac
np.matmul(jac.T, jac)
```

or, perhaps faster but python 3 only:

```python
hess = ret.jac.T @ ret.jac
if scipy.sparse.issparse(hess):
    hess = hess.todense()
```

#### A Minimal, Complete, and Verifiable example

```python
import lmfit

def f(params):
    return params['x'] - 1, params['y'] - 1

params = lmfit.Parameters()
params.add('x', value=0)
params.add('y', value=0)
out = lmfit.minimize(f, params, method='least_squares', jac_sparsity=[[1, 1], [1, 1]])
```

#### Error message:

```
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/Admin/.pyenv/versions/3.7.3/Python.framework/Versions/3.7/lib/python3.7/site-packages/lmfit/minimizer.py", line 2448
, in minimize
    return fitter.minimize(method=method)
  File "/Users/Admin/.pyenv/versions/3.7.3/Python.framework/Versions/3.7/lib/python3.7/site-packages/lmfit/minimizer.py", line 2117
, in minimize
    return function(**kwargs)
  File "/Users/Admin/.pyenv/versions/3.7.3/Python.framework/Versions/3.7/lib/python3.7/site-packages/lmfit/minimizer.py", line 1430
, in least_squares
    hess = np.matmul(ret.jac.T, ret.jac)
ValueError: matmul: Input operand 0 does not have enough dimensions (has 0, gufunc core with signature (n?,k),(k,m?)->(n?,m?) requi
res 1)
```

#### Version information

```
Python: 3.7.3 (default, May 26 2019, 09:41:28)
[Clang 10.0.1 (clang-1001.0.46.4)]
lmfit: 0.9.14, scipy: 1.3.0, numpy: 1.16.3, asteval: 0.9.14, uncertainties: 3.1, six: 1.12.0
```
diff --git a/tests/test_least_squares.py b/tests/test_least_squares.py
index 778c2049..b9aa64e0 100644
--- a/tests/test_least_squares.py
+++ b/tests/test_least_squares.py
@@ -2,6 +2,8 @@
 import numpy as np
 from numpy.testing import assert_allclose
 import pytest
+from scipy.sparse import bsr_matrix
+from scipy.sparse.linalg import aslinearoperator
 
 import lmfit
 from lmfit.models import VoigtModel
@@ -116,3 +118,52 @@ def test_least_squares_solver_options(peakdata, capsys):
 
     assert 'Iteration' in captured.out
     assert 'final cost' in captured.out
+
+
+def test_least_squares_jacobian_types():
+    """Test support for Jacobian of all types supported by least_squares."""
+    # Build function
+    # f(x, y) = (x - a)^2 + (y - b)^2
+    np.random.seed(42)
+    a = np.random.normal(0, 1, 50)
+    np.random.seed(43)
+    b = np.random.normal(0, 1, 50)
+
+    def f(params):
+        return (params['x'] - a)**2 + (params['y'] - b)**2
+
+    # Build analytic Jacobian functions with the different possible return types
+    # 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))
+
+    def jac_sparse(params, *args, **kwargs):
+        return bsr_matrix(jac_array(params, *args, **kwargs))
+
+    def jac_operator(params, *args, **kwargs):
+        return aslinearoperator(jac_array(params, *args, **kwargs))
+    # Build parameters
+    params = lmfit.Parameters()
+    params.add('x', value=0)
+    params.add('y', value=0)
+    # Solve model for numerical Jacobian and each analytic Jacobian function
+    result = lmfit.minimize(f, params, method='least_squares')
+    result_array = lmfit.minimize(
+        f, params, method='least_squares',
+        jac=jac_array)
+    result_sparse = lmfit.minimize(
+        f, params, method='least_squares',
+        jac=jac_sparse)
+    result_operator = lmfit.minimize(
+        f, params, method='least_squares',
+        jac=jac_operator)
+    # Check that all have uncertainties
+    assert result.errorbars
+    assert result_array.errorbars
+    assert result_sparse.errorbars
+    assert result_operator.errorbars
+    # Check that all have ~equal covariance matrix
+    assert_allclose(result.covar, result_array.covar)
+    assert_allclose(result.covar, result_sparse.covar)
+    assert_allclose(result.covar, result_operator.covar)

pip install -e .[dev]
parse_log_pytest
requirements.txt
pytest coverage codecov
apt-get update apt-get install -y gcc
3.7
requirements.txt
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
asteval==0.9.31
certifi @ file:///croot/certifi_1671487769961/work/certifi
charset-normalizer==3.4.1
codecov==2.1.13
coverage==7.2.7
exceptiongroup==1.2.2
future==1.0.0
idna==3.10
importlib-metadata==6.7.0
iniconfig==2.0.0
-e git+https://github.com/lmfit/lmfit-py.git@ffe550c9c1e3f2c6b82b7175dd101f5c6cd8f685#egg=lmfit
numpy==1.21.6
packaging==24.0
pluggy==1.2.0
pytest==7.4.4
requests==2.31.0
scipy==1.7.3
six==1.17.0
tomli==2.0.1
typing_extensions==4.7.1
uncertainties==3.1.7
urllib3==2.0.7
zipp==3.15.0

name: lmfit-py
channels:
  - defaults
  - https://repo.anaconda.com/pkgs/main
  - https://repo.anaconda.com/pkgs/r
  - conda-forge
dependencies:
  - _libgcc_mutex=0.1=main
  - _openmp_mutex=5.1=1_gnu
  - ca-certificates=2025.2.25=h06a4308_0
  - certifi=2022.12.7=py37h06a4308_0
  - ld_impl_linux-64=2.40=h12ee557_0
  - libffi=3.4.4=h6a678d5_1
  - libgcc-ng=11.2.0=h1234567_1
  - libgomp=11.2.0=h1234567_1
  - libstdcxx-ng=11.2.0=h1234567_1
  - ncurses=6.4=h6a678d5_0
  - openssl=1.1.1w=h7f8727e_0
  - pip=22.3.1=py37h06a4308_0
  - python=3.7.16=h7a1cb2a_0
  - readline=8.2=h5eee18b_0
  - setuptools=65.6.3=py37h06a4308_0
  - sqlite=3.45.3=h5eee18b_0
  - tk=8.6.14=h39e8969_0
  - wheel=0.38.4=py37h06a4308_0
  - xz=5.6.4=h5eee18b_1
  - zlib=1.2.13=h5eee18b_1
  - pip:
      - asteval==0.9.31
      - charset-normalizer==3.4.1
      - codecov==2.1.13
      - coverage==7.2.7
      - exceptiongroup==1.2.2
      - future==1.0.0
      - idna==3.10
      - importlib-metadata==6.7.0
      - iniconfig==2.0.0
      - numpy==1.21.6
      - packaging==24.0
      - pluggy==1.2.0
      - pytest==7.4.4
      - requests==2.31.0
      - scipy==1.7.3
      - six==1.17.0
      - tomli==2.0.1
      - typing-extensions==4.7.1
      - uncertainties==3.1.7
      - urllib3==2.0.7
      - zipp==3.15.0
prefix: /opt/conda/envs/lmfit-py

tests/test_least_squares.py::test_least_squares_jacobian_types
tests/test_least_squares.py::test_least_squares_with_bounds tests/test_least_squares.py::test_least_squares_cov_x[False] tests/test_least_squares.py::test_least_squares_cov_x[True] tests/test_least_squares.py::test_least_squares_solver_options
lmfit/minimizer.py
lmfit/minimizer.py