ENH: handle empty matrices in qr decomposition #11593
Conversation
numpy/linalg/linalg.py
Outdated
| if _isEmpty2d(a): | ||
| k = min(m, n) | ||
| if mode == 'reduced': | ||
| # ‘reduced’ : returns q, r with dimensions (M, K), (K, N) (default) |
There was a problem hiding this comment.
The forward quotes are unicode, causing a test failure
There was a problem hiding this comment.
oh... actually intended to remove those lines; leaving them there and fixing the quotes.
|
Needs at least a passing mention in the docs, and in the release notes. It would be cool to add the gufunc signatures for all the linalg routines to the docstrings too, but that is probably a different issue. |
|
My sense would be that supporting empty matrices should be the expectation, and that non-support should be what is flagged. It seems that this is consistent with |
|
But that said, what should be documented and where? |
|
|
|
Done. |
doc/release/1.16.0-notes.rst
Outdated
| Previously, a ``LinAlgError`` would be raised when empty matrix ("flat" | ||
| or "skinny") is passed in. This has been fixed so that outputs of | ||
| appropriate shapes are returned for the various modes. | ||
|
|
There was a problem hiding this comment.
I wonder if this and the similar note about lstsq for your related work in #11594 could be combined since they are both about handling empty matrices for linalg functions, but maybe that's a pain to coordinate between 2 separate PRs & could be done after I suppose.
There was a problem hiding this comment.
Separated them because for lstsq, I was working in a strange way (site-packages). For qr it was a more straightforward "change and run tests".
|
What do you mean in two places? |
|
Once in the size == 0 path, and once in the regular path. The only place that 0 needs to be a special case is when actually filling the arrays - the existing code already makes all the arrays the right shape, it just sometimes forgets to fill them, leaving them uninitialized |
|
Looking further into it, I also note that the But I'm not inclined to change that (even though the tests would pass) because it might break something and benefit nothing (given deprecation). |
…case of empty arrays (they do not play well with empty arrays)
Codecov Report
@@ Coverage Diff @@
## master #11593 +/- ##
=========================================
Coverage ? 85.7%
=========================================
Files ? 327
Lines ? 82002
Branches ? 0
=========================================
Hits ? 70281
Misses ? 11721
Partials ? 0
Continue to review full report at Codecov.
|
|
Patch looks pretty good now, thanks - I do wonder if we can just call the lapack functions though - often the only issue is that we didn't read the docs carefully on how to pass size-0 parameters |
|
Let me take another swing at that. I just happened to be looking at the LAPACK documentation earlier today and I'm going have another look. I don't like "not calling anything". |
numpy/linalg/linalg.py
Outdated
| if results['info'] != 0: | ||
| raise LinAlgError('%s returns %d' % (routine_name, results['info'])) | ||
|
|
||
| # compute q | ||
| lwork = int(abs(work[0])) | ||
| lwork = max(1, int(abs(work[0]))) |
There was a problem hiding this comment.
I'd be surprised if this is necessary for the work spaces, since lapack should be populating them with whatever it really needs
There was a problem hiding this comment.
It's not. But my problem was having to pass in lwork=0 which was invalid (it had to be at least... oh shit I may have made a mistake...)
http://www.netlib.org/lapack/explore-3.1.1-html/dgeqrf.f.html
http://www.netlib.org/lapack/explore-3.1.1-html/zungqr.f.html
Tests pass but I'm going to go fix it...
| @@ -875,14 +875,14 @@ def qr(a, mode='reduced'): | |||
| # calculate optimal size of work data 'work' | |||
| lwork = 1 | |||
| work = zeros((lwork,), t) | |||
| results = lapack_routine(m, n, a, m, tau, work, -1, 0) | |||
| results = lapack_routine(m, n, a, max(1, m), tau, work, -1, 0) | |||
There was a problem hiding this comment.
This is exactly the fix I was expecting :)
There was a problem hiding this comment.
... yes... and (the late) Gene Golub is not going to tell me to give back my e-mail account. (I got my gmail invite from him.)
numpy/linalg/tests/test_linalg.py
Outdated
| assert_almost_equal(np.dot(q, r), a) | ||
| assert_almost_equal(np.dot(q.T.conj(), q), np.eye(k)) | ||
| assert_almost_equal(np.triu(r), r) | ||
| # |
There was a problem hiding this comment.
No need for empty comments - just leave a blank line
numpy/linalg/tests/test_linalg.py
Outdated
| @@ -1583,8 +1583,45 @@ def check_qr(self, a): | |||
| assert_almost_equal(r2, r1) | |||
|
|
|||
| def test_qr_empty(self): | |||
| a = np.zeros((0, 2)) | |||
| assert_raises(linalg.LinAlgError, linalg.qr, a) | |||
| for m, n in [(3, 0), (0, 3)]: | |||
There was a problem hiding this comment.
Can you add (0, 0) here too?
numpy/linalg/tests/test_linalg.py
Outdated
| assert_(r.shape == (k, n)) | ||
| assert_almost_equal(np.dot(q, r), a) | ||
| assert_almost_equal(np.dot(q.T.conj(), q), np.eye(k)) | ||
| assert_almost_equal(np.triu(r), r) |
There was a problem hiding this comment.
I'm not sure these three tests are interesting - in these tests k == 0, a.size == 0, and r.size == 0 - so there are no values to compare anyway
There was a problem hiding this comment.
I'll drop some of the ones that are implied
There was a problem hiding this comment.
Actually - just call self.check_qr(a), and it will do all the work for you
|
See my comment above - almost all of your test can be replaced with a call to |
numpy/linalg/tests/test_linalg.py
Outdated
|
|
||
| self.check_qr(a) | ||
|
|
||
| r = np.linalg.qr(a, mode='r') |
There was a problem hiding this comment.
Isn't this already covered by check_qr?
|
Great. Non-ugly patches are good. |
|
@convexset: I think you might have missed my one lingering comment above - one of your tests seems to still duplicate |
|
@eric-wieser: it's not exactly duplicating |
|
Can you elaborate on in what way it's not duplicated? Does your test test more than what check_qr tests? If so, can you move your extra checks to within |
numpy/linalg/tests/test_linalg.py
Outdated
| r = np.linalg.qr(a, mode='r') | ||
| assert_equal(r.dtype, a_dtype) | ||
| assert_(isinstance(r, a_type)) | ||
| assert_equal(r.shape, (k, n)) |
There was a problem hiding this comment.
This shape test should probably be in check_qr anyway
There was a problem hiding this comment.
It's there implicitly.
# mode == 'reduced'
q1, r1 = linalg.qr(a, mode='reduced')
# ...
assert_(r1.shape == (k, n))
# ...
# mode == 'r'
r2 = linalg.qr(a, mode='r')
# ...
assert_almost_equal(r2, r1)So I took that bit out.
There was a problem hiding this comment.
Fair enough. I think the explicit shape check you had here was a little clearer, but it's not super important.
numpy/linalg/tests/test_linalg.py
Outdated
| assert_(isinstance(r, a_type)) | ||
| assert_equal(r.shape, (k, n)) | ||
|
|
||
| h, tau = np.linalg.qr(a, mode='raw') |
There was a problem hiding this comment.
This test can stay - for whatever reason, it seems check_qr decided not to test it
numpy/linalg/tests/test_linalg.py
Outdated
| (4, 0, 1), | ||
| (4, 0, 2), | ||
| (4, 2, 0), | ||
| (0, 0, 0) |
There was a problem hiding this comment.
Theses are triplets, but you only have two variables. I assume you just want [(4, 0), (0, 4), (0, 0)]?
|
This comment still applies |
|
Thanks for your persistence and your first contribution, @convexset ! |
handle empty matrices in qr decomposition; added tests