# BUG: "sp.linalg.solve_discrete_are" fails for random data #6572

Merged
merged 20 commits into from Sep 19, 2016
+352 −85

## Conversation

Member

### ilayn commented Sep 13, 2016

 This fixes #2251 and its tests per the discussion on the issue page. The doc string is updated.
 BUG: "sp.linalg.solve_discrete_are" fails for random data 
 72bd91b 
This fixes #2251 and its tests per the discussion on the issue page. The
doc string is updated.
mentioned this pull request Sep 13, 2016
added 6 commits Sep 13, 2016
 fixed the unpacking shortcut 
 fd231f7 
This is Py2 issue I guess
 Fixed the wrong number of arguments. 
 d3feed6 
Hopefully this will make Travis happy.
 made the test residual explicit 
 afd16ab 
 branch conflict resolved (again) 
 5120f0c 
 switching to arrays 
 90fc11d 
For some unknown reason, the test works locally but trips up on the
master branch. Here testing whether this is due to the matrix object.
 Finally found the embarassing omission 
 d3c3f56 
Member

### person142 commented Sep 13, 2016

 BTW you can check all the pep8 failures locally by installing pep8 (pip install pep8) and running pep8 scipy from the top-level scipy directory. It's looking like that will be the last thing you need for the Travis build to succeed.
added 3 commits Sep 13, 2016
 PEP8 clean-up 
 f9162f8 
 Revert "PEP8 clean-up" 
 922e57e 
This reverts commit f9162f8.
 PEP8 clean-up 
 97eb01b 
Member Author

### ilayn commented Sep 13, 2016

 Bah, I guess I'm not cut for this. I've checked the files thrice with pep8 and it still fails.
 PEP8 cleanup also for the tests 
 4058607 
Member Author

### ilayn commented Sep 14, 2016

 I've slept over it checked with some more examples and I think this PR is ready for review. Thanks .
 U.S. Energy Research and Development Agency under contract ERDA-E(49-18)-2087. http://dspace.mit.edu/bitstream/handle/1721.1/1301/R-0859-05666488.pdf P. van Dooren , "A Generalized Eigenvalue Approach For Solving Riccati

#### person142 Sep 14, 2016

Member

I know it was like this in the original, but it would be better to put this in a references section and just cite it here. See:

https://github.com/numpy/numpy/blob/master/doc/HOWTO_DOCUMENT.rst.txt#docstring-standard

 @@ -304,12 +310,20 @@ def solve_discrete_are(a, b, q, r): .. math:: X = A'XA-(A'XB)(R+B'XB)^-1(B'XA)+Q It is solved directly using a Schur decomposition method. It is solved via forming the extended symplectic pencil

#### person142 Sep 14, 2016

Member

Probably would fit better in the notes section.

 .. math:H - \lambda J given by the block matrices [ A 0 B ] [ I 0 B ]

#### person142 Sep 14, 2016

Member

Should this just define H = [...], J = [...] instead of restating the equation? I guess you could also take out the first H - \lambda J.

#### ilayn Sep 14, 2016

Author Member

That's a domain specific formulation. I've moved all of those into Notes.

 if n != r.shape[0]: raise ValueError("Matrix b and r should have the same number of cols.") # Check if the data is (sufficiently) hermitian

#### person142 Sep 14, 2016

Member

I don't think other linalg routines do checks like this, so I wouldn't put them here. But I might be wrong.

#### ilayn Sep 14, 2016

Author Member

If we don't do it, it will trip somewhere in the code with the cryptic broadcast exception. I think it is better to show why it fails such that the user immediately spots the mistake.

 # Shape consistency checks m, n = b.shape

#### person142 Sep 14, 2016

Member

Nitpick: no extra line here or between the ifs.

#### person142 Sep 18, 2016

Member

Don't forget the ifs.

 overwrite_b=True, check_finite=False, output=out_str) # Get the relevant parts of the stable subspace basis # TODO: Check if it succeded to get the right ordering.

#### person142 Sep 14, 2016

Member

Could you say more about the TODO? Are there places where this could fail?

#### ilayn Sep 14, 2016

Author Member

These equations don't have solvability guarantees. Depending on the data matrices a solution might not exist. I don't know what SciPy way when no solution is found : return empty matrix or raise exception ?

 np.array([[-2j,], [1j,]]), np.array([[1, 0], [0, 2]]), np.array([[1,],])), # Real a, b; complex q, r (corrected for hermitian Q)

#### person142 Sep 14, 2016

Member

Safe to take out the "corrected for hermitian Q" bit--if someone can't see this diff they won't know what it means.

Member

### person142 commented Sep 14, 2016

 Though the previous method doesn't always work, it still might be nice to see how the two compare speed-wise. (Especially on more sizable problems.) Benchmarks would go here: https://github.com/scipy/scipy/blob/master/benchmarks/benchmarks/linalg.py You can run them with python runtests.py --bench (you will need asv; pip install asv). An intro to writing them is here: https://asv.readthedocs.io/en/latest/writing_benchmarks.html
Member Author

### ilayn commented Sep 14, 2016

 Unfortunately, on a windows machine I can't do any of the build related stuff including the runtests. But since the previous one was not reliable a speedtest won't make too much of a difference. There are many extras to this implementation to make it comparable to matlab. But this version is at least consistent.
Member

### person142 commented Sep 15, 2016 • edited

 But since the previous one was not reliable a speedtest won't make too much of a difference. I guess I worry that someone will have working code that suddenly slows down a lot. Breaking their case to make other cases work might not make them happy. If there's a significant speed difference it's possible that we might want to have both versions provided that for the old one there's a way to either document a reasonable subset of problems on which it will work or detect failures and raise an exception. Maybe I'm just being silly though. Edit: yes, I think I'm being silly.
Member

### person142 commented Sep 15, 2016

 Also it's possible to run benchmarks without runtests.py by doing something like asv dev --bench  which won't trigger a build and only runs the benchmark once. (So take the results with a grain of salt.)
Member

### person142 commented Sep 15, 2016

 It would be nice to include more DAREX examples in the tests.
Member Author

### ilayn commented Sep 15, 2016

 Indeed but then I have to include balancing algorithms which is basically copying everything from my other implementation which doesn't seem that easy to convert to SciPy way without breaking some back-compatibility.
Member Author

### ilayn commented Sep 17, 2016 • edited

 @person142 Oh that I can tell. In the CAREX examples, there are two examples fail** even after my balancing implementations (ex.6 and 12 in the matlab version). In DAREX I have to finish the conversion from matlab to Python but so far so good but I'm sure it will fail on a few suspicious looking ones. One major pain is the PEP8 and those matrices don't go well together, is there a way to ignore PEP8 temporarily of things or what is the preferred way of entering them? Tall Submatrices ? ** By fail I mean the accuracy is not good enough absolutely though satisfactory in the relative sense (max error is less than a certain percent of the 1-norm the matrix).
Member

### pv commented Sep 17, 2016 • edited

 If the matrices are big, I would put them as data files (e.g. text files), under scipy/linalg/tests/data, and load via mat = np.loadtxt(os.path.join(os.path.dirname(__file__), 'data', 'darex_1.dat')) etc. (or, alternatively, put as npy data files, anything goes) If they are small, I would just follow pep8 format --- reformatting probably doesn't take that much time even if done manually.
 added DAREX library examples to tests 
 24f1250 
removed previous two test cases as they are already in darex, fails are
commented out
Member Author

### ilayn commented Sep 17, 2016

 darex #12 is what I mean by having a relatively accurate but not acceptable in terms of absolute terms solution, as the solution should be diag(1,1e12) but gets approx. diag(1,1e12 + 2.5e7).
Member

### person142 commented Sep 18, 2016

 Ok, for the benefit of anyone looking at this: with the new set of tests the old function has 4 exceptions due to ill-conditioning in r and one failure. See e.g. here: https://github.com/person142/scipy/tree/old-discrete-are
requested changes
Member

### person142 left a comment

 Ok, I left some comments. They're mostly pretty small things; I think this is pretty close to being done.
 ERDA-E(49-18)-2087. http://dspace.mit.edu/bitstream/handle/1721.1/1301/R-0859-05666488.pdf The equation is solved by forming the extended symplectic matrix pencil, as described in [1]_, .. math:H - \lambda J given by the block matrices

#### person142 Sep 18, 2016

Member
 [ A 0 B ] [ I 0 B ] [ -Q I 0 ] - \lambda * [ 0 A^T 0 ] [ 0 0 R ] [ 0 -B^T 0 ]

#### person142 Sep 18, 2016

Member

This isn't rendering right.

#### pv Sep 18, 2016

Member

Needs :: (two colons) at the end of the preceding text, ie. ... given by the block matrices::

 if np.iscomplexobj(x): r_or_c = complex if not np.equal(*x.shape):

#### person142 Sep 18, 2016

Member

This fails for a 1 x 1 array.

 if m != a.shape[0]: raise ValueError("Matrix a and b should have the same number of rows.")

Member

#### ilayn Sep 18, 2016

Author Member

I didn't get your point on spacing

#### person142 Sep 18, 2016

Member

It's not a show-stopper, but I'd like related groups of if statements to not have any spaces between them (and maybe be converted to if/elif's) so that they cluster together visually. But this might just be a personal preference so I won't press the point.

#### ilayn Sep 18, 2016

Author Member

I thought you wanted a spacing between the ifs inside the for loop so I extrapolated from that, but I have no preference.

 if m != q.shape[0]: raise ValueError("Matrix a and q should have the same shape.")

#### person142 Sep 18, 2016

Member

Ditto on spacing.

 up, ul, uu = lu(u00) if 1/cond(uu) < np.spacing(1.): return np.array([[]])

#### person142 Sep 18, 2016 • edited

Member

Sorry, forgot to clarify this before: failure conditions should raise an exception (in this case a LinAlgError).

 sym_threshold = np.max(np.spacing([1000., norm(u_sym, 1)])) if np.max(np.abs(u_sym)) > sym_threshold: return np.array([[]])

#### person142 Sep 18, 2016

Member

Exception here too.

 [3.7730, -30.2800, 14.6900]]) * 0.001, np.diag([50, 0, 0, 0, 50, 0, 0, 0, 0]), np.eye(3)), ## darex #12

#### person142 Sep 18, 2016

Member

I have some suggested changes for this test:

https://github.com/person142/scipy/blob/ilayn/scipy/linalg/tests/test_solvers.py#L136

It turns the failing tests into known failures so that they're more visible.

reviewed
 np.array([[1, 2], [1, 3]]), np.array([[1, 1+1j], [1-1j, 2]]), np.array([[2, -2j], [2j, 3]])), # User-reported

#### person142 Sep 18, 2016

Member

For got to mention this: reference the original github issue here so it's easy to find later.

#### ilayn Sep 18, 2016

Author Member

Ah, forgot to include this...

added 4 commits Sep 18, 2016
 now raises exceptions and 2d array 
 77a3edb 
also cleanup and fixes in the doc string
 Revert "now raises exceptions and 2d array" 
 419257e 
This reverts commit 77a3edb.
 Revert "Revert "now raises exceptions and 2d array"" 
 05fcb4c 
This reverts commit 419257e.
 minor pep8 
 3ddfa5f 
reviewed
 np.eye(3), 1e6 * np.eye(3), 1e6 * np.eye(3), "Fails"),

#### person142 Sep 18, 2016

Member

Ok so I got a little lazy putting in the causes for failure here; should probably at least be something like "unsatisfactory absolute error" if that's the cause.

#### ilayn Sep 18, 2016

Author Member

No actually it is a complete failure to find a valid solution. The other two I marked it as bad accuracy

 explicit fail reasons and github issue number 
 918cf45 
Member Author

### ilayn commented Sep 19, 2016

 This hopefully concludes the requested changes. Please let me know if there are additional items
approved these changes
Member

### person142 commented Sep 19, 2016

 I'm happy (the CircleCI passed, not sure why it's still yellow), but I'll wait to get a second opinion from another developer.
reviewed
 # Check the deviation from symmetry for success u_sym = u00.conj().T.dot(u10) u_sym -= u_sym.conj().T

#### pv Sep 19, 2016

Member

This will bomb if u_sym is real, because in-place operations on overlapping memory regions is undefined behavior.
For extra fun, you don't see it for small matrices due to buffering.

#### pv Sep 19, 2016

Member

Iow, should be u_sym = u_sym - u_sym.conj().T

#### ilayn Sep 19, 2016

Author Member

Wow. That's really good to know. I'll fix it in a moment. Is there any place where I can read more about it?

#### pv Sep 19, 2016 • edited

Member

It's the issue about Numpy arrays and views to numpy arrays, eg b = a[::2]; b[0] = 1 modifies also a. I think it's mentioned in most numpy tutorials. With in-place operations etc., you want to make sure the output operand is not a view of any kind of the input.

However, I think this particular caveat here will go away in near future numpy/numpy#8043. There are also long threads about this in the numpy issue tracker, but I'd recommend those only for entertainment purposes now that there's a consensus on how the issue will be addressed.

Member

### pv commented Sep 19, 2016

 Looks OK to me apart from that one cmment.
Member

### person142 commented Sep 19, 2016

 Also would be good to squash things.
 combined ifs, changed loop var x to a geniune name 
 b56d107 
blank line reduction, removal of redundant indexing
Member

### pv commented Sep 19, 2016

 @person142: there's an option to squash and merge in the github merge button
Member Author

### ilayn commented Sep 19, 2016

 I don't know what "squash" refers to actually.
merged commit 9eaeddc into scipy:master Sep 19, 2016
3 checks passed
3 checks passed
ci/circleci Your tests passed on CircleCI!
Details
codecov/project 78.79% (target 1.00%)
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
Member

### pv commented Sep 19, 2016

 Thanks, merged!
added this to the 0.19.0 milestone Sep 19, 2016
Member

### person142 commented Sep 19, 2016

 Great! Thanks for sticking with us @ilayn, and I look forward to your future improvements!
Member Author

### ilayn commented Sep 19, 2016

 Pleasure is mine.
This was referenced Sep 23, 2016