LSMR fails on some random values when providing an x0

[gabrielbriones]

In my computer the isolve.test() fails in testing lsmr() at this point: https://github.com/scipy/scipy/blob/master/scipy/sparse/linalg/isolve/tests/test_lsmr.py#L85 due to x_ref and x arrays are not almost equal.

Since the matrix G and array b are randomly defined in lsqr.py I decided to test the code running a stand alone script multiple times, assuming that it will generate a different G and b values each time.

$ cat test_lsmr.py
from scipy.sparse.linalg.isolve.tests.test_lsqr import G, b
from scipy.sparse.linalg import lsmr
from numpy.testing import assert_array_almost_equal

x_ref = lsmr(G, b)[0]
x0 = lsmr(G, b, maxiter=1)[0]
x = lsmr(G, b, x0=x0)[0]
assert_array_almost_equal(x_ref, x)

I can see that sometimes the test code passes and sometimes it did not:

$ python3 test_lsmr.py
$ python3 test_lsmr.py
$ python3 test_lsmr.py
$ python3 test_lsmr.py
$ python3 test_lsmr.py
$ python3 test_lsmr.py
$ python3 test_lsmr.py
$ python3 test_lsmr.py
$ python3 test_lsmr.py
$ python3 test_lsmr.py
Traceback (most recent call last):
  File "test_lsmr.py", line 8, in <module>
    assert_array_almost_equal(x_ref, x)
  File "/usr/lib/python3.6/site-packages/numpy/testing/utils.py", line 962, in assert_array_almost_equal
    precision=decimal)
  File "/usr/lib/python3.6/site-packages/numpy/testing/utils.py", line 778, in assert_array_compare
    raise AssertionError(msg)
AssertionError: 
Arrays are not almost equal to 6 decimals

(mismatch 77.14285714285714%)
 x: array([ 1.343296,  0.604053, -0.111929, -0.125148,  0.714433,  1.80346 ,
       -0.618127, -0.146294, -1.349041,  1.810285, -1.841975,  0.004248,
        0.077046,  0.720121,  0.523681,  0.074032,  1.277396, -0.143822,...
 y: array([ 1.343293,  0.604054, -0.111933, -0.125147,  0.714436,  1.803458,
       -0.618128, -0.146297, -1.349044,  1.810282, -1.841972,  0.004249,
        0.077043,  0.720116,  0.523683,  0.074031,  1.277398, -0.143826,...

$ python3 test_lsmr.py
$ python3 test_lsmr.py
Traceback (most recent call last):
  File "test_lsmr.py", line 8, in <module>
    assert_array_almost_equal(x_ref, x)
  File "/usr/lib/python3.6/site-packages/numpy/testing/utils.py", line 962, in assert_array_almost_equal
    precision=decimal)
  File "/usr/lib/python3.6/site-packages/numpy/testing/utils.py", line 778, in assert_array_compare
    raise AssertionError(msg)
AssertionError: 
Arrays are not almost equal to 6 decimals

(mismatch 11.42857142857143%)
 x: array([ 1.466697, -1.815673, -0.059471, -1.066072,  1.518105,  0.40146 ,
        0.57259 ,  0.245936,  1.166329, -2.271805, -0.087285,  0.483411,
        0.659616,  1.291392,  1.082163,  0.324475, -1.165986,  0.719474,...
 y: array([ 1.466698, -1.815674, -0.059473, -1.066073,  1.518104,  0.40146 ,
        0.572589,  0.245936,  1.166328, -2.271805, -0.087286,  0.483411,
        0.659616,  1.291391,  1.082163,  0.324475, -1.165986,  0.719472,...

I also tried modifying the size of G and b, I realize that it is reproducible only when size >= 25 (The default is 35) and becomes more often in sizes bigger than 55. I saved the values of G and b when the test fails and hardcoded it into the script.
My assumption is that the default tolerance is not appropiate.

Reproducing code example:

from scipy.sparse.linalg import lsmr
from numpy.testing import assert_array_almost_equal
from numpy import array

G = array([[4.5176785744636865, 11.701926270820824, 7.913875419887336, 24.689044975219378, 4.055144379698328, 1.003473942762554, 9.087566188989971, 6.670303889071847, 4.872305855451147, 5.127624634910172, 7.037884654464798, 12.959097577599596, 1.570440577219601, 8.894190270248613, 3.2325095814238196, 1.2036534035056083, 3.8490798680699, 3.5814271086917664, 3.798808372546076, 8.48210758147563, 5.990069833761698, 2.879403611202188, 3.305935490636988, 3.055946392829132, 1.7878289129306917], [3.517678574463688, 12.701926270820824, 7.913875419887336, 24.689044975219378, 4.055144379698328, 1.003473942762554, 9.087566188989971, 6.670303889071847, 4.872305855451147, 5.127624634910172, 7.037884654464798, 12.959097577599596, 1.570440577219601, 8.894190270248613, 3.2325095814238196, 1.2036534035056083, 3.8490798680699, 3.5814271086917664, 3.798808372546076, 8.48210758147563, 5.990069833761698, 2.879403611202188, 3.305935490636988, 3.055946392829132, 1.7878289129306917], [3.517678574463688, 11.701926270820824, 8.913875419887336, 24.689044975219378, 4.055144379698328, 1.003473942762554, 9.087566188989971, 6.670303889071847, 4.872305855451147, 5.127624634910172, 7.037884654464798, 12.959097577599596, 1.570440577219601, 8.894190270248613, 3.2325095814238196, 1.2036534035056083, 3.8490798680699, 3.5814271086917664, 3.798808372546076, 8.48210758147563, 5.990069833761698, 2.879403611202188, 3.305935490636988, 3.055946392829132, 1.7878289129306917], [3.517678574463688, 11.701926270820824, 7.913875419887336, 25.689044975219378, 4.055144379698328, 1.003473942762554, 9.087566188989971, 6.670303889071847, 4.872305855451147, 5.127624634910172, 7.037884654464798, 12.959097577599596, 1.570440577219601, 8.894190270248613, 3.2325095814238196, 1.2036534035056083, 3.8490798680699, 3.5814271086917664, 3.798808372546076, 8.48210758147563, 5.990069833761698, 2.879403611202188, 3.305935490636988, 3.055946392829132, 1.7878289129306917], [3.517678574463688, 11.701926270820824, 7.913875419887336, 24.689044975219378, 5.055144379698328, 1.003473942762554, 9.087566188989971, 6.670303889071847, 4.872305855451147, 5.127624634910172, 7.037884654464798, 12.959097577599596, 1.570440577219601, 8.894190270248613, 3.2325095814238196, 1.2036534035056083, 3.8490798680699, 3.5814271086917664, 3.798808372546076, 8.48210758147563, 5.990069833761698, 2.879403611202188, 3.305935490636988, 3.055946392829132, 1.7878289129306917], [3.517678574463688, 11.701926270820824, 7.913875419887336, 24.689044975219378, 4.055144379698328, 2.0034739427625547, 9.087566188989971, 6.670303889071847, 4.872305855451147, 5.127624634910172, 7.037884654464798, 12.959097577599596, 1.570440577219601, 8.894190270248613, 3.2325095814238196, 1.2036534035056083, 3.8490798680699, 3.5814271086917664, 3.798808372546076, 8.48210758147563, 5.990069833761698, 2.879403611202188, 3.305935490636988, 3.055946392829132, 1.7878289129306917], [3.517678574463688, 11.701926270820824, 7.913875419887336, 24.689044975219378, 4.055144379698328, 1.003473942762554, 10.087566188989971, 6.670303889071847, 4.872305855451147, 5.127624634910172, 7.037884654464798, 12.959097577599596, 1.570440577219601, 8.894190270248613, 3.2325095814238196, 1.2036534035056083, 3.8490798680699, 3.5814271086917664, 3.798808372546076, 8.48210758147563, 5.990069833761698, 2.879403611202188, 3.305935490636988, 3.055946392829132, 1.7878289129306917], [3.517678574463688, 11.701926270820824, 7.913875419887336, 24.689044975219378, 4.055144379698328, 1.003473942762554, 9.087566188989971, 7.670303889071847, 4.872305855451147, 5.127624634910172, 7.037884654464798, 12.959097577599596, 1.570440577219601, 8.894190270248613, 3.2325095814238196, 1.2036534035056083, 3.8490798680699, 3.5814271086917664, 3.798808372546076, 8.48210758147563, 5.990069833761698, 2.879403611202188, 3.305935490636988, 3.055946392829132, 1.7878289129306917], [3.517678574463688, 11.701926270820824, 7.913875419887336, 24.689044975219378, 4.055144379698328, 1.003473942762554, 9.087566188989971, 6.670303889071847, 5.872305855451147, 5.127624634910172, 7.037884654464798, 12.959097577599596, 1.570440577219601, 8.894190270248613, 3.2325095814238196, 1.2036534035056083, 3.8490798680699, 3.5814271086917664, 3.798808372546076, 8.48210758147563, 5.990069833761698, 2.879403611202188, 3.305935490636988, 3.055946392829132, 1.7878289129306917], [3.517678574463688, 11.701926270820824, 7.913875419887336, 24.689044975219378, 4.055144379698328, 1.003473942762554, 9.087566188989971, 6.670303889071847, 4.872305855451147, 6.127624634910172, 7.037884654464798, 12.959097577599596, 1.570440577219601, 8.894190270248613, 3.2325095814238196, 1.2036534035056083, 3.8490798680699, 3.5814271086917664, 3.798808372546076, 8.48210758147563, 5.990069833761698, 2.879403611202188, 3.305935490636988, 3.055946392829132, 1.7878289129306917], [3.517678574463688, 11.701926270820824, 7.913875419887336, 24.689044975219378, 4.055144379698328, 1.003473942762554, 9.087566188989971, 6.670303889071847, 4.872305855451147, 5.127624634910172, 8.037884654464797, 12.959097577599596, 1.570440577219601, 8.894190270248613, 3.2325095814238196, 1.2036534035056083, 3.8490798680699, 3.5814271086917664, 3.798808372546076, 8.48210758147563, 5.990069833761698, 2.879403611202188, 3.305935490636988, 3.055946392829132, 1.7878289129306917], [3.517678574463688, 11.701926270820824, 7.913875419887336, 24.689044975219378, 4.055144379698328, 1.003473942762554, 9.087566188989971, 6.670303889071847, 4.872305855451147, 5.127624634910172, 7.037884654464798, 13.959097577599598, 1.570440577219601, 8.894190270248613, 3.2325095814238196, 1.2036534035056083, 3.8490798680699, 3.5814271086917664, 3.798808372546076, 8.48210758147563, 5.990069833761698, 2.879403611202188, 3.305935490636988, 3.055946392829132, 1.7878289129306917], [3.517678574463688, 11.701926270820824, 7.913875419887336, 24.689044975219378, 4.055144379698328, 1.003473942762554, 9.087566188989971, 6.670303889071847, 4.872305855451147, 5.127624634910172, 7.037884654464798, 12.959097577599596, 2.5704405772196006, 8.894190270248613, 3.2325095814238196, 1.2036534035056083, 3.8490798680699, 3.5814271086917664, 3.798808372546076, 8.48210758147563, 5.990069833761698, 2.879403611202188, 3.305935490636988, 3.055946392829132, 1.7878289129306917], [3.517678574463688, 11.701926270820824, 7.913875419887336, 24.689044975219378, 4.055144379698328, 1.003473942762554, 9.087566188989971, 6.670303889071847, 4.872305855451147, 5.127624634910172, 7.037884654464798, 12.959097577599596, 1.570440577219601, 9.894190270248615, 3.2325095814238196, 1.2036534035056083, 3.8490798680699, 3.5814271086917664, 3.798808372546076, 8.48210758147563, 5.990069833761698, 2.879403611202188, 3.305935490636988, 3.055946392829132, 1.7878289129306917], [3.517678574463688, 11.701926270820824, 7.913875419887336, 24.689044975219378, 4.055144379698328, 1.003473942762554, 9.087566188989971, 6.670303889071847, 4.872305855451147, 5.127624634910172, 7.037884654464798, 12.959097577599596, 1.570440577219601, 8.894190270248613, 4.232509581423819, 1.2036534035056083, 3.8490798680699, 3.5814271086917664, 3.798808372546076, 8.48210758147563, 5.990069833761698, 2.879403611202188, 3.305935490636988, 3.055946392829132, 1.7878289129306917], [3.517678574463688, 11.701926270820824, 7.913875419887336, 24.689044975219378, 4.055144379698328, 1.003473942762554, 9.087566188989971, 6.670303889071847, 4.872305855451147, 5.127624634910172, 7.037884654464798, 12.959097577599596, 1.570440577219601, 8.894190270248613, 3.2325095814238196, 2.2036534035056086, 3.8490798680699, 3.5814271086917664, 3.798808372546076, 8.48210758147563, 5.990069833761698, 2.879403611202188, 3.305935490636988, 3.055946392829132, 1.7878289129306917], [3.517678574463688, 11.701926270820824, 7.913875419887336, 24.689044975219378, 4.055144379698328, 1.003473942762554, 9.087566188989971, 6.670303889071847, 4.872305855451147, 5.127624634910172, 7.037884654464798, 12.959097577599596, 1.570440577219601, 8.894190270248613, 3.2325095814238196, 1.2036534035056083, 4.8490798680699, 3.5814271086917664, 3.798808372546076, 8.48210758147563, 5.990069833761698, 2.879403611202188, 3.305935490636988, 3.055946392829132, 1.7878289129306917], [3.517678574463688, 11.701926270820824, 7.913875419887336, 24.689044975219378, 4.055144379698328, 1.003473942762554, 9.087566188989971, 6.670303889071847, 4.872305855451147, 5.127624634910172, 7.037884654464798, 12.959097577599596, 1.570440577219601, 8.894190270248613, 3.2325095814238196, 1.2036534035056083, 3.8490798680699, 4.5814271086917655, 3.798808372546076, 8.48210758147563, 5.990069833761698, 2.879403611202188, 3.305935490636988, 3.055946392829132, 1.7878289129306917], [3.517678574463688, 11.701926270820824, 7.913875419887336, 24.689044975219378, 4.055144379698328, 1.003473942762554, 9.087566188989971, 6.670303889071847, 4.872305855451147, 5.127624634910172, 7.037884654464798, 12.959097577599596, 1.570440577219601, 8.894190270248613, 3.2325095814238196, 1.2036534035056083, 3.8490798680699, 3.5814271086917664, 4.798808372546076, 8.48210758147563, 5.990069833761698, 2.879403611202188, 3.305935490636988, 3.055946392829132, 1.7878289129306917], [3.517678574463688, 11.701926270820824, 7.913875419887336, 24.689044975219378, 4.055144379698328, 1.003473942762554, 9.087566188989971, 6.670303889071847, 4.872305855451147, 5.127624634910172, 7.037884654464798, 12.959097577599596, 1.570440577219601, 8.894190270248613, 3.2325095814238196, 1.2036534035056083, 3.8490798680699, 3.5814271086917664, 3.798808372546076, 9.48210758147563, 5.990069833761698, 2.879403611202188, 3.305935490636988, 3.055946392829132, 1.7878289129306917], [3.517678574463688, 11.701926270820824, 7.913875419887336, 24.689044975219378, 4.055144379698328, 1.003473942762554, 9.087566188989971, 6.670303889071847, 4.872305855451147, 5.127624634910172, 7.037884654464798, 12.959097577599596, 1.570440577219601, 8.894190270248613, 3.2325095814238196, 1.2036534035056083, 3.8490798680699, 3.5814271086917664, 3.798808372546076, 8.48210758147563, 6.990069833761699, 2.879403611202188, 3.305935490636988, 3.055946392829132, 1.7878289129306917], [3.517678574463688, 11.701926270820824, 7.913875419887336, 24.689044975219378, 4.055144379698328, 1.003473942762554, 9.087566188989971, 6.670303889071847, 4.872305855451147, 5.127624634910172, 7.037884654464798, 12.959097577599596, 1.570440577219601, 8.894190270248613, 3.2325095814238196, 1.2036534035056083, 3.8490798680699, 3.5814271086917664, 3.798808372546076, 8.48210758147563, 5.990069833761698, 3.8794036112021875, 3.305935490636988, 3.055946392829132, 1.7878289129306917], [3.517678574463688, 11.701926270820824, 7.913875419887336, 24.689044975219378, 4.055144379698328, 1.003473942762554, 9.087566188989971, 6.670303889071847, 4.872305855451147, 5.127624634910172, 7.037884654464798, 12.959097577599596, 1.570440577219601, 8.894190270248613, 3.2325095814238196, 1.2036534035056083, 3.8490798680699, 3.5814271086917664, 3.798808372546076, 8.48210758147563, 5.990069833761698, 2.879403611202188, 4.305935490636989, 3.055946392829132, 1.7878289129306917], [3.517678574463688, 11.701926270820824, 7.913875419887336, 24.689044975219378, 4.055144379698328, 1.003473942762554, 9.087566188989971, 6.670303889071847, 4.872305855451147, 5.127624634910172, 7.037884654464798, 12.959097577599596, 1.570440577219601, 8.894190270248613, 3.2325095814238196, 1.2036534035056083, 3.8490798680699, 3.5814271086917664, 3.798808372546076, 8.48210758147563, 5.990069833761698, 2.879403611202188, 3.305935490636988, 4.0559463928291315, 1.7878289129306917], [3.517678574463688, 11.701926270820824, 7.913875419887336, 24.689044975219378, 4.055144379698328, 1.003473942762554, 9.087566188989971, 6.670303889071847, 4.872305855451147, 5.127624634910172, 7.037884654464798, 12.959097577599596, 1.570440577219601, 8.894190270248613, 3.2325095814238196, 1.2036534035056083, 3.8490798680699, 3.5814271086917664, 3.798808372546076, 8.48210758147563, 5.990069833761698, 2.879403611202188, 3.305935490636988, 3.055946392829132, 2.7878289129306917]])
b = array([0.014327751645286208, 0.3705206141618929, -0.19738512561554636, -0.25040840806933995, -1.1802300643746382, 0.7684470150116808, 0.7109157743111, 2.3331717375390397, 1.3868109880876647, 1.0861940449952845, 0.18820958017195968, 1.2059588092447053, 0.7043244743077872, 0.8085336500332501, 1.8874926911143695, 0.5621026091687727, -1.324228685719128, -0.8510221053727485, 0.8027201320897752, -1.0995904838280939, -1.4488597214826695, 0.35394687611316356, -0.5897253716602514, 1.3167421451497685, -1.3047822231094595])

x_ref = lsmr(G, b)[0]
x0 = lsmr(G, b, maxiter=1)[0]
x = lsmr(G, b, x0=x0)[0]
assert_array_almost_equal(x_ref, x)

Error message:

Traceback (most recent call last):
  File "fail_lsmr.py", line 11, in <module>
    assert_array_almost_equal(x_ref, x)
  File "/usr/lib/python3.6/site-packages/numpy/testing/utils.py", line 962, in assert_array_almost_equal
    precision=decimal)
  File "/usr/lib/python3.6/site-packages/numpy/testing/utils.py", line 778, in assert_array_compare
    raise AssertionError(msg)
AssertionError: 
Arrays are not almost equal to 6 decimals

(mismatch 4.0%)
 x: array([-0.228055,  0.128139, -0.439767, -0.492789, -1.422613,  0.526064,
        0.468534,  2.090789,  1.144428,  0.843812, -0.054173,  0.963577,
        0.461942,  0.566152,  1.64511 ,  0.31972 , -1.566611, -1.093405,...
 y: array([-0.228054,  0.128139, -0.439767, -0.49279 , -1.422612,  0.526065,
        0.468534,  2.09079 ,  1.144429,  0.843812, -0.054172,  0.963577,
        0.461943,  0.566152,  1.645111,  0.319721, -1.566611, -1.093404,...

Scipy/Numpy/Python version information:

1.0.0 1.13.3 sys.version_info(major=3, minor=6, micro=3, releaselevel='final', serial=0)

[VictorRodriguez]

docker run -it clearlinux/machine-learning
gets you a local docker (with shell into it) that has all the bits needed to repro

[ilayn]

Could you also add which architecture are you running the tests, 32- or 64- bit, which Linux and so on?

[pv]

assert_almost_equal tests 1e-7 precision but the default precision for lsmr is 1e-6, so I think it is not expected the solutions are the same.

In particular, it's well known that floating-point reproducibility even for the same problem is not obtained with modern optimized linear algebra libraries.

[pv]

I guess the test should either set atol/btol appropriately or relax the assert_close tolerance

[gabrielbriones]

It's a 64-bit Clearlinux 19090 with kernel 4.13 and Intel Core i5-6260U. I have also seem the same behavior in Fedora 26, but the failing values of G and b across distributions. I had to run the test code multiple times in Fedora to get failing values of G and b.

[ilayn]

We should also fix the random seed to be able to reproduce such problems easily.