ENH: handle empty matrices in qr decomposition #11593

convexset · 2018-07-20T09:46:56Z

handle empty matrices in qr decomposition; added tests

mattip · 2018-07-20T14:32:40Z

numpy/linalg/linalg.py

+    if _isEmpty2d(a):
+        k = min(m, n)
+        if mode == 'reduced':
+            # ‘reduced’ : returns q, r with dimensions (M, K), (K, N) (default)


The forward quotes are unicode, causing a test failure

oh... actually intended to remove those lines; leaving them there and fixing the quotes.

mattip · 2018-07-20T14:39:19Z

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.

convexset · 2018-07-20T14:41:17Z

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 dot, tensordot and other operations.

convexset · 2018-07-20T14:41:52Z

But that said, what should be documented and where?

mattip · 2018-07-20T16:11:28Z

doc/release/1.16.0-notes.rst under Changes or Improvements.

convexset · 2018-07-20T16:17:02Z

Done.

tylerjereddy · 2018-07-20T16:47:54Z

doc/release/1.16.0-notes.rst

+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.
+


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.

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".

eric-wieser

We shouldn't allocate the output arrays in two places - the arrays allocated in the normal case are already the right size and shape. It should be sufficient to add q[...] = eye, to correct for anything lapack fails to do.

convexset · 2018-07-22T09:49:21Z

What do you mean in two places?

eric-wieser · 2018-07-22T20:50:03Z

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

convexset · 2018-07-29T06:29:52Z

Looking further into it, I also note that the raw and economic modes don't match up with the documentation.

    if mode == 'raw':
        return a, tau

    if mode == 'economic':
        if t != result_t :
            a = a.astype(result_t, copy=False)
        return wrap(a.T)

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-io · 2018-07-29T07:13:34Z

Codecov Report

❗ No coverage uploaded for pull request base (master@977431a). Click here to learn what that means.
The diff coverage is 100%.

@@            Coverage Diff            @@
##             master   #11593   +/-   ##
=========================================
  Coverage          ?    85.7%           
=========================================
  Files             ?      327           
  Lines             ?    82002           
  Branches          ?        0           
=========================================
  Hits              ?    70281           
  Misses            ?    11721           
  Partials          ?        0

Impacted Files	Coverage Δ
numpy/linalg/linalg.py	`91.05% <100%> (ø)`
numpy/linalg/tests/test_linalg.py	`97.11% <100%> (ø)`

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 977431a...ff9063c. Read the comment docs.

eric-wieser · 2018-07-29T16:39:12Z

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

New approach taken looks good

convexset · 2018-07-29T17:38:15Z

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".

eric-wieser · 2018-07-29T17:50:34Z

numpy/linalg/linalg.py

    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])))


I'd be surprised if this is necessary for the work spaces, since lapack should be populating them with whatever it really needs

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...

eric-wieser · 2018-07-29T17:50:52Z

numpy/linalg/linalg.py

@@ -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)


This is exactly the fix I was expecting :)

... 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.)

eric-wieser

I suspect you'll still need a case to populate the identity matrix when the input is empty but the output is not

eric-wieser · 2018-07-29T18:43:21Z

numpy/linalg/tests/test_linalg.py

+            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)
+            #


No need for empty comments - just leave a blank line

eric-wieser · 2018-07-29T18:44:04Z

numpy/linalg/tests/test_linalg.py

@@ -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)]:


Can you add (0, 0) here too?

eric-wieser · 2018-07-29T18:45:19Z

numpy/linalg/tests/test_linalg.py

+            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)


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

I'll drop some of the ones that are implied

Actually - just call self.check_qr(a), and it will do all the work for you

eric-wieser · 2018-07-29T18:55:15Z

See my comment above - almost all of your test can be replaced with a call to self.check_qr(a)

eric-wieser · 2018-07-29T19:19:14Z

numpy/linalg/tests/test_linalg.py

+
+            self.check_qr(a)
+
+            r = np.linalg.qr(a, mode='r')


Isn't this already covered by check_qr?

convexset · 2018-07-29T20:03:05Z

Great. Non-ugly patches are good.

eric-wieser · 2018-07-31T06:27:24Z

@convexset: I think you might have missed my one lingering comment above - one of your tests seems to still duplicate check_qr - or am I missing something?

convexset · 2018-07-31T06:49:46Z

@eric-wieser: it's not exactly duplicating check_qr... but it is implied...

eric-wieser · 2018-07-31T06:51:01Z

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 check_qr?

eric-wieser · 2018-07-31T06:52:12Z

numpy/linalg/tests/test_linalg.py

+            r = np.linalg.qr(a, mode='r')
+            assert_equal(r.dtype, a_dtype)
+            assert_(isinstance(r, a_type))
+            assert_equal(r.shape, (k, n))


This shape test should probably be in check_qr anyway

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.

Fair enough. I think the explicit shape check you had here was a little clearer, but it's not super important.

eric-wieser · 2018-07-31T06:52:40Z

numpy/linalg/tests/test_linalg.py

+            assert_(isinstance(r, a_type))
+            assert_equal(r.shape, (k, n))
+
+            h, tau = np.linalg.qr(a, mode='raw')


This test can stay - for whatever reason, it seems check_qr decided not to test it

eric-wieser · 2018-07-31T06:54:30Z

numpy/linalg/tests/test_linalg.py

+        (4, 0, 1),
+        (4, 0, 2),
+        (4, 2, 0),
+        (0, 0, 0)


Theses are triplets, but you only have two variables. I assume you just want [(4, 0), (0, 4), (0, 0)]?

eric-wieser · 2018-07-31T06:55:10Z

This comment still applies

eric-wieser

Commits need squashing, but I can do that when I merge. Let's hope tests pass!

eric-wieser · 2018-07-31T17:00:35Z

Thanks for your persistence and your first contribution, @convexset !

handle empty matrices in qr decomposition

36fcb0c

This was referenced Jul 20, 2018

np.linalg.qr: handle empty matrices #11412

Closed

ENH: Implement np.linalg.{svd,qr,lstsq} for matrices with .size == 0 #8654

Closed

mattip reviewed Jul 20, 2018

View reviewed changes

fixed offending quotes

3fa693a

mattip added 01 - Enhancement component: numpy.linalg 56 - Needs Release Note. Needs an entry in doc/release/upcoming_changes labels Jul 20, 2018

release notes updated

6dabc39

tylerjereddy reviewed Jul 20, 2018

View reviewed changes

tylerjereddy changed the title ~~handle empty matrices in qr decomposition~~ ENH: handle empty matrices in qr decomposition Jul 20, 2018

mattip approved these changes Jul 22, 2018

View reviewed changes

eric-wieser previously requested changes Jul 22, 2018

View reviewed changes

revised to use existing array sizes; skipping calls to lapack in the …

448b409

…case of empty arrays (they do not play well with empty arrays)

setting things back where they were and calling LAPACK functions better

1624a80

eric-wieser reviewed Jul 29, 2018

View reviewed changes

updated tests

7f22736

simplified tests

68d608e

eric-wieser reviewed Jul 29, 2018

View reviewed changes

eric-wieser removed the 56 - Needs Release Note. Needs an entry in doc/release/upcoming_changes label Jul 31, 2018

eric-wieser added this to the 1.16.0 release milestone Jul 31, 2018

eric-wieser reviewed Jul 31, 2018

View reviewed changes

removed redundant test; expressed sub-test with pytest.mark.parametrize

5809e9a

eric-wieser reviewed Jul 31, 2018

View reviewed changes

fixed dumb copy-and-paste error

16d79e9

eric-wieser approved these changes Jul 31, 2018

View reviewed changes

convexset added 2 commits July 31, 2018 15:18

clarified terminology

c4c6426

trimmed excess on line

ff9063c

eric-wieser merged commit 8fdc446 into numpy:master Jul 31, 2018

convexset deleted the handle-edge-cases branch August 1, 2018 05:41

mattip mentioned this pull request Aug 7, 2018

WIP: Allow lstsq on empty arrays #11604

Closed

Uh oh!

ENH: handle empty matrices in qr decomposition #11593

ENH: handle empty matrices in qr decomposition #11593

Uh oh!

Conversation

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Codecov Report

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!