Repository: scipy/scipy
Pull request: BUG: fix small tolerances handling for minpack and add a test.
#### Reference issue
fix #11572 

#### What does this implement/fix?
fix small tolerances handling for minpack and add a test.

#### Additional information
The root cause of this issue is that, minimum tolerance of minpack is bigger than computer EPS.
In my environment, EPS is 2e-16, but minpack cannot handle tolerances which is smaller than 1e-10.
So, I added a new minpack status for too small tolerances, which is based on:
Maxima Manual: 61. minpack http://maxima.sourceforge.net/docs/manual/de/maxima_61.html

Problem statement:
least_squares: too small tolerances not catched with method='lm'
<!-- 
Thank you for taking the time to file a bug report. 
Please fill in the fields below, deleting the sections that 
don't apply to your issue. You can view the final output
by clicking the preview button above.
Note: This is a comment, and won't appear in the output.
-->

If I pass very small but not too small termination tolerances to `least_squares` with `method='lm'`, instead of raising an exception immediately, the routine is called, it catches the error and `least_squares` does not handle it properly.

#### Reproducing code example:
<!-- 
If you place your code between the triple backticks below, 
it will be rendered as a code block. 
-->

```python
from scipy import optimize
import numpy as np

def try_fit(y, H, p0, **kwargs):
    def res(p):
        return y - H @ p
    try:
        result = optimize.least_squares(res, p0, ftol=0, gtol=0, **kwargs)
    except ValueError:
        pass # "at least one of the tolerances must be > 2e-16"

p0 = np.array([-0.62631426,  0.34736521, -1.53904034])
y = np.array([ 0.21042216, -0.774209  , -1.69254418, -0.20248775,  2.15496739,
               0.57563308,  2.14036788,  1.9187491 , -1.03636067,  3.58819572])
H = np.array([[ 0.4620464 , -1.91682697, -0.55389804],
              [-0.79682045, -1.60787256,  1.58394725],
              [-0.52023376, -1.67694436,  0.18293149],
              [ 0.59203689, -1.37410441, -0.34793584],
              [-0.76610329,  0.39296886, -0.79426522],
              [-0.0211633 ,  1.20145896,  0.22471291],
              [-0.09212487,  0.35775539, -0.083304  ],
              [-0.91542477,  0.1543938 , -1.14452311],
              [ 0.53072674,  1.28778141,  0.73426207],
              [-0.71480276,  0.66167675, -0.94704367]])

try_fit(y, H, p0, xtol=0, method='lm')
try_fit(y, H, p0, xtol=1e-13, method='trf')
try_fit(y, H, p0, xtol=1e-16, method='lm')
try_fit(y, H, p0, xtol=1e-13, method='lm')
```

#### Error message:
<!-- If any, paste the *full* error message inside a code block
as above (starting from line Traceback)
-->

```python
KeyError                                  Traceback (most recent call last)
~/Documents/prog/yalsi/scipybug.py in <module>
     27 try_fit(y, H, p0, xtol=1e-13, method='trf')
     28 try_fit(y, H, p0, xtol=1e-16, method='lm')
---> 29 try_fit(y, H, p0, xtol=1e-13, method='lm')

~/Documents/prog/yalsi/scipybug.py in try_fit(y, H, p0, **kwargs)
      6         return y - H @ p
      7     try:
----> 8         result = optimize.least_squares(res, p0, ftol=0, gtol=0, **kwargs)
      9     except ValueError:
     10         pass # "at least one of the tolerances must be > 2e-16"

~/Documents/Programmi/python3_virtualenv/lib/python3.6/site-packages/scipy/optimize/_lsq/least_squares.py in least_squares(fun, x0, jac, bounds, method, ftol, xtol, gtol, x_scale, loss, f_scale, diff_step, tr_solver, tr_options, jac_sparsity, max_nfev, verbose, args, kwargs)
    910     if method == 'lm':
    911         result = call_minpack(fun_wrapped, x0, jac_wrapped, ftol, xtol, gtol,
--> 912                               max_nfev, x_scale, diff_step)
    913 
    914     elif method == 'trf':

~/Documents/Programmi/python3_virtualenv/lib/python3.6/site-packages/scipy/optimize/_lsq/least_squares.py in call_minpack(fun, x0, jac, ftol, xtol, gtol, max_nfev, x_scale, diff_step)
     87     njev = info.get('njev', None)
     88 
---> 89     status = FROM_MINPACK_TO_COMMON[status]
     90     active_mask = np.zeros_like(x0, dtype=int)
     91 

KeyError: 6
```

#### Scipy/Numpy/Python version information:
```
1.4.1 1.16.4 sys.version_info(major=3, minor=6, micro=4, releaselevel='final', serial=0)
```

Hints from issue discussion:
@Gattocrucco 

What if you try to fix it? It looks like a good easy issue to become familiar with scipy development process.

Commit message: fix tolerance check for method lm.

Patch to review:
diff --git a/scipy/optimize/_lsq/least_squares.py b/scipy/optimize/_lsq/least_squares.py
index c70656c0cde7..9e8ab211742e 100644
--- a/scipy/optimize/_lsq/least_squares.py
+++ b/scipy/optimize/_lsq/least_squares.py
@@ -102,7 +102,7 @@ def prepare_bounds(bounds, n):
     return lb, ub
 
 
-def check_tolerance(ftol, xtol, gtol):
+def check_tolerance(ftol, xtol, gtol, method):
     def check(tol, name):
         if tol is None:
             tol = 0
@@ -116,7 +116,10 @@ def check(tol, name):
     xtol = check(xtol, "xtol")
     gtol = check(gtol, "gtol")
 
-    if ftol < EPS and xtol < EPS and gtol < EPS:
+    if method == "lm" and (ftol < EPS or xtol < EPS or gtol < EPS):
+        raise ValueError("All tolerances must be higher than machine epsilon"
+                         "({:.2e}) for method 'lm'.".format(EPS))
+    elif ftol < EPS and xtol < EPS and gtol < EPS:
         raise ValueError("At least one of the tolerances must be higher than "
                          "machine epsilon ({:.2e}).".format(EPS))
 
@@ -793,7 +796,7 @@ def least_squares(
 
     x_scale = check_x_scale(x_scale, x0)
 
-    ftol, xtol, gtol = check_tolerance(ftol, xtol, gtol)
+    ftol, xtol, gtol = check_tolerance(ftol, xtol, gtol, method)
 
     def fun_wrapped(x):
         return np.atleast_1d(fun(x, *args, **kwargs))
diff --git a/scipy/optimize/tests/test_least_squares.py b/scipy/optimize/tests/test_least_squares.py
index 041439419c06..0d1bdff409c0 100644
--- a/scipy/optimize/tests/test_least_squares.py
+++ b/scipy/optimize/tests/test_least_squares.py
@@ -4,7 +4,7 @@
 from numpy.linalg import norm
 from numpy.testing import (assert_, assert_allclose,
                            assert_equal, suppress_warnings)
-from pytest import raises as assert_raises
+from pytest import raises as assert_raises, warns
 from scipy.sparse import issparse, lil_matrix
 from scipy.sparse.linalg import aslinearoperator
 
@@ -369,6 +369,8 @@ def test_error_raised_when_all_tolerances_below_eps(self):
                       method=self.method, ftol=None, xtol=None, gtol=None)
 
     def test_convergence_with_only_one_tolerance_enabled(self):
+        if self.method == 'lm':
+            return # should not do test
         x0 = [-2, 1]
         x_opt = [1, 1]
         for ftol, xtol, gtol in [(1e-8, None, None),
@@ -744,3 +746,30 @@ def test_basic():
     # test that 'method' arg is really optional
     res = least_squares(fun_trivial, 2.0)
     assert_allclose(res.x, 0, atol=1e-10)
+
+
+def test_small_tolerances_for_lm():
+    def res(p):
+        return y - H @ p
+
+    p0 = np.array([-0.62631426, 0.34736521, -1.53904034])
+    y = np.array([0.21042216, -0.774209, -1.69254418, -0.20248775, 2.15496739,
+                  0.57563308, 2.14036788, 1.9187491, -1.03636067, 3.58819572])
+    H = np.array([[0.4620464, -1.91682697, -0.55389804],
+                  [-0.79682045, -1.60787256, 1.58394725],
+                  [-0.52023376, -1.67694436, 0.18293149],
+                  [0.59203689, -1.37410441, -0.34793584],
+                  [-0.76610329, 0.39296886, -0.79426522],
+                  [-0.0211633, 1.20145896, 0.22471291],
+                  [-0.09212487, 0.35775539, -0.083304],
+                  [-0.91542477, 0.1543938, -1.14452311],
+                  [0.53072674, 1.28778141, 0.73426207],
+                  [-0.71480276, 0.66167675, -0.94704367]])
+
+    for ftol, xtol, gtol in [(0.0, 1e-13, 1e-13),
+                             (1e-13, 0.0, 1e-13),
+                             (1e-13, 1e-13, 0.0)]:
+        with warns(UserWarning): # UseWarning for EPS are shown.
+            print(ftol, xtol, gtol)
+            assert_raises(ValueError, least_squares, res, p0, xtol=xtol,
+                      ftol=ftol, gtol=gtol, method='lm')


Code review comments:
diff_hunk: @@ -744,3 +746,30 @@ def test_basic():
     # test that 'method' arg is really optional
     res = least_squares(fun_trivial, 2.0)
     assert_allclose(res.x, 0, atol=1e-10)
+
+
+def test_small_tolerances_for_lm():
+    def res(p):
+        return y - H @ p
+
+    p0 = np.array([-0.62631426, 0.34736521, -1.53904034])
+    y = np.array([0.21042216, -0.774209, -1.69254418, -0.20248775, 2.15496739,
+                  0.57563308, 2.14036788, 1.9187491, -1.03636067, 3.58819572])
+    H = np.array([[0.4620464, -1.91682697, -0.55389804],
+                  [-0.79682045, -1.60787256, 1.58394725],
+                  [-0.52023376, -1.67694436, 0.18293149],
+                  [0.59203689, -1.37410441, -0.34793584],
+                  [-0.76610329, 0.39296886, -0.79426522],
+                  [-0.0211633, 1.20145896, 0.22471291],
+                  [-0.09212487, 0.35775539, -0.083304],
+                  [-0.91542477, 0.1543938, -1.14452311],
+                  [0.53072674, 1.28778141, 0.73426207],
+                  [-0.71480276, 0.66167675, -0.94704367]])
+
+    for ftol, xtol, gtol in [(0.0, 1e-13, 1e-13),
path: scipy/optimize/tests/test_least_squares.py
text: @user1:
I suggest that you should use `None` instead of `0.0` here and remove `warns` context.

@author:
I fix it.

diff_hunk: @@ -744,3 +746,30 @@ def test_basic():
     # test that 'method' arg is really optional
     res = least_squares(fun_trivial, 2.0)
     assert_allclose(res.x, 0, atol=1e-10)
+
+
+def test_small_tolerances_for_lm():
+    def res(p):
+        return y - H @ p
+
+    p0 = np.array([-0.62631426, 0.34736521, -1.53904034])
+    y = np.array([0.21042216, -0.774209, -1.69254418, -0.20248775, 2.15496739,
+                  0.57563308, 2.14036788, 1.9187491, -1.03636067, 3.58819572])
+    H = np.array([[0.4620464, -1.91682697, -0.55389804],
+                  [-0.79682045, -1.60787256, 1.58394725],
+                  [-0.52023376, -1.67694436, 0.18293149],
+                  [0.59203689, -1.37410441, -0.34793584],
+                  [-0.76610329, 0.39296886, -0.79426522],
+                  [-0.0211633, 1.20145896, 0.22471291],
+                  [-0.09212487, 0.35775539, -0.083304],
+                  [-0.91542477, 0.1543938, -1.14452311],
+                  [0.53072674, 1.28778141, 0.73426207],
+                  [-0.71480276, 0.66167675, -0.94704367]])
+
+    for ftol, xtol, gtol in [(0.0, 1e-13, 1e-13),
+                             (1e-13, 0.0, 1e-13),
+                             (1e-13, 1e-13, 0.0)]:
+        with warns(UserWarning): # UseWarning for EPS are shown.
+            print(ftol, xtol, gtol)
path: scipy/optimize/tests/test_least_squares.py
text: @user1:
Debug print.

@author:
I removed it.