Repository: scipy/scipy
Pull request: Enable setting minres initial vector
#### Reference issue
Closes https://github.com/scipy/scipy/issues/6843

#### What does this implement/fix?
Makes sure `scipy.sparse.linalg.isolve.minres` uses parameter `x0` correctly.

#### Additional information
Combines https://github.com/scipy/scipy/pull/6887 and https://github.com/scipy/scipy/pull/7080

Tests for `x0` are added to `test_minres` not to `test_iterative` where they can be tested for more solvers. A first attempt gave some errors and I'm not familiar enough with all solvers to fix it immediately. If requested, I can look into it.

Problem statement:
minres always uses zero vector as initial vector and user doesn't know about it
What I mean, is _scipy/sparse/linalg/isolve/minres_, which is in python just `scipy.sparse.linalg.minres`.

In 46 line:
`    A,M,x,b,postprocess = make_system(A,M,x0,b,xtype)`

It is usual preprocessing, which just check is `x0` equals to `None`. But in 88 line we have:
`    x = zeros(n, dtype=xtype)`

So, anyway initial vector is zero. Is there a reason to allow user to specify `x0` ? It doesn't make a sense. Maybe better delete `x0=None` from the function signature?

Hints from issue discussion:
Maybe better to do it like this:

```
def minres(A, b, shift=0.0, tol=1e-5, maxiter=None, xtype=None,
           M=None, callback=None, show=False, check=False):
    A,M,x,b,postprocess = make_system(A,M,None,b,xtype)
```


Is there a reason why you want to eliminate `x0` instead of making it work as advertised?
It doesn't work as advertised: it doesn't uses `x0` anyway. All code just rewrited directly from corresponding Fortran files (you can find them here: http://web.stanford.edu/group/SOL/software/minres/), where were no need in initial vector (cause they assign it to zero vector).

In my opinion the problem is that nobody says that it is. There isn't any mentions that `x0` is useless in that particular function. Original signature just copypasted from other functions. Because of that some stupid students like me tries to use this function and they wonder why it works as it works, not as they expect (and this expectations are provided by wrong signature).
> It doesn't work as advertised: it doesn't uses `x0` anyway

Right, we're agreed on that. The question is whether `x0` should be eliminated or whether the code should be rewritten to actually use `x0`.
i.e. one could deprecate `x0` in the signature, or one could remove the `x = zeros(n, dtype=xtype)` line. (As long is there isn't a mathematical reason for not doing that.)
I don't see the reason to deprecate x0. Presumably, minres can work with given initial vector also.
Yeah, that's my feeling too.
I just noticed this issue.  There is no reason why x0 should or should not be 0. 
The MinRes method converges for all initial guess x0, if A is (real) positive definite.

But the documentation tells us:  "A can be indefinite or singular."

Furthermore, if I can read your code correctly you even check if A is symmetric:
[here](https://github.com/scipy/scipy/blob/master/scipy/sparse/linalg/isolve/minres.py#L111) and [there](https://github.com/scipy/scipy/blob/master/scipy/sparse/linalg/isolve/minres.py#L171)

Which is not related, since
A=[1 , 1 ; -1 , 1]  
is positive definite (since x^t A x = x₀² + x₁² > 0, ∀x=(x₀,x₁) =/= 0 ) ,
but not symmetric.

Oh and about line 88: x = zeros(...)
We have in line 94:  r1 = b
This actually should be  r1 = b-Ax, if you intend to use the input.  
"x=zeros" is used nowhere explcitly, only implicitly in line 94.
I can confirm this bug.

Commit message: activated test_gh_6843, added new test for correct x0

Patch to review:
diff --git a/scipy/sparse/linalg/isolve/minres.py b/scipy/sparse/linalg/isolve/minres.py
index 22f6808c9a46..846a3102fb74 100644
--- a/scipy/sparse/linalg/isolve/minres.py
+++ b/scipy/sparse/linalg/isolve/minres.py
@@ -110,23 +110,19 @@ def minres(A, b, x0=None, shift=0.0, tol=1e-5, maxiter=None,
 
     eps = finfo(xtype).eps
 
-    x = zeros(n, dtype=xtype)
-
     # Set up y and v for the first Lanczos vector v1.
     # y  =  beta1 P' v1,  where  P = C**(-1).
     # v is really P' v1.
 
-    y = b
-    r1 = b
-
-    y = psolve(b)
+    r1 = b - A*x
+    y = psolve(r1)
 
-    beta1 = inner(b,y)
+    beta1 = inner(y, y)
 
     if beta1 < 0:
         raise ValueError('indefinite preconditioner')
     elif beta1 == 0:
-        return (postprocess(x), 0)
+        return (postprocess(zeros(n, dtype=xtype)), 0)
 
     beta1 = sqrt(beta1)
 
diff --git a/scipy/sparse/linalg/isolve/tests/test_minres.py b/scipy/sparse/linalg/isolve/tests/test_minres.py
index b94304bb0d8a..c0e30215f2d2 100644
--- a/scipy/sparse/linalg/isolve/tests/test_minres.py
+++ b/scipy/sparse/linalg/isolve/tests/test_minres.py
@@ -3,7 +3,8 @@
 import numpy as np
 from numpy.testing import assert_equal, assert_allclose, assert_
 from scipy.sparse.linalg.isolve import minres
-import pytest
+from scipy.linalg import norm
+
 from pytest import raises as assert_raises
 from .test_iterative import assert_normclose
 
@@ -27,9 +28,7 @@ def test_singular():
     assert_normclose(A.dot(xp), b, tol=1e-5)
 
 
-@pytest.mark.skip(reason="Skip Until gh #6843 is fixed")
-def test_gh_6843():
-    """check if x0 is being used by tracing iterates"""
+def test_x0_is_used_by():
     A, b = get_sample_problem()
     # Random x0 to feed minres
     np.random.seed(12345)
@@ -63,3 +62,14 @@ def test_asymmetric_fail():
     A[2, 1] = 2
     with assert_raises(ValueError):
         xp, info = minres(A, b, check=True)
+
+
+def test_minres_non_default_x0():
+    np.random.seed(1234)
+    tol = 10**(-6)
+    a = np.random.randn(5, 5)
+    a = np.dot(a, a.T)
+    b = np.random.randn(5)
+    c = np.random.randn(5)
+    x = minres(a, b, x0=c, tol=tol)[0]
+    assert norm(a.dot(x) - b) < tol


Code review comments:
diff_hunk: @@ -110,23 +110,19 @@ def minres(A, b, x0=None, shift=0.0, tol=1e-5, maxiter=None,
 
     eps = finfo(xtype).eps
 
-    x = zeros(n, dtype=xtype)
-
     # Set up y and v for the first Lanczos vector v1.
     # y  =  beta1 P' v1,  where  P = C**(-1).
     # v is really P' v1.
 
-    y = b
-    r1 = b
-
-    y = psolve(b)
+    r1 = b - A*x
+    y = psolve(r1)
 
-    beta1 = inner(b,y)
+    beta1 = inner(y, y)
path: scipy/sparse/linalg/isolve/minres.py
text: @user1:
Did you check that this change is correct? It's now different from previous (with preconditioning), also for `x0=0`.

@author:
I was indeed to focused on no preconditioning, let me get back to it and start by adding more tests!

@author:
Corrected to a version that works also with preconditioning and added a test

@user1:
And then a few more questions: do you understand the minres algorithm? Did you check that these changes are correct (without assuming the previous PRs are)?

@author:
I have quite some experience and understanding of (lanczos bases) optimisation algorithms. However, in practice, all is limited to non-preconditioned systems!

In a Lanczos based approach, if you have an initial vector, basically all you need to change is the initial residual, which happens on line 117 in the PR. This is also compatible with the way the original code deals with next residuals: https://github.com/scipy/scipy/blob/4d26a150340a51153b01615c8afbe433380340d8/scipy/sparse/linalg/isolve/minres.py#L196
Notice that the code now also reproduces exactly what happened before in case `x0=0`.
I looked in the literature and as far as I can tell, this introduction of `x0` is indeed correct. However, as I said, I have no extended practical experience with preconditioners and given that I'm working on someone else code, I cannot say I'm 100% sure that there is subtle mistake, sorry. 

@user1:
Thanks, was just asking before thinking it through myself. After reading it properly, the changes look ok.