Add routines for updating QR decompositions #4249
Conversation
|
I don't have much idea about the details. One question: The update routines require In the stats applications Another question: Is there anything about singular matrices? How does it behave? Is it supposed to work or is it ruled out? Is it stable in the near singular case, or is it better to regularize in that case on the user side? |
|
And, Thanks for working on this. There will be many use cases where this is useful to speed up the calculations for example in statsmodels. |
|
|
||
| p_subdiag_qr(m, n-p, q, qs, r, rs, k, p, work) | ||
|
|
||
| libc.stdlib.free(work) |
There was a problem hiding this comment.
Cython docs, http://docs.cython.org/src/tutorial/memory_allocation.html, seem to recommend a try-finally block for free. For my education: is it that it's actually not needed?
There was a problem hiding this comment.
I'm not actually sure. It seems to me that if it does, there must be tons of wrapped code that can leak memory in certain circumstances (multiple threads etc.), given that this compiles down to a simple C function. Perhaps some one that knows for certain can weigh in?
There was a problem hiding this comment.
Depends on whether something in between can raise or not. In this case,
since it's a nogil function, no exceptions can be raised.
|
Nothing works with economy decompositions right now. I don't know of any reason why it couldn't be made to work though. I stopped here mostly because this is already 2000+ lines and that seemed big enough for a first cut at this. There is nothing here that does anything in particular with singular matrices. Do you have anything in particular in mind? Or a reference? |
Ok, understandable. The increased memory handling might eat away all advantages for efficient updatting for most of our (statsmodels) application.
This was mainly a usage question. 2 issues when I was reading around (skimming many articles) for this. |
|
I've pushed support for rank-k updates to economic mode decompositions. Right now, the update is |
|
Can these be generalized to handle |
|
I suppose it could be. In that case I'd end up scaling the rows of The general approach to performing updates is to reduce |
|
Yup, you're right of course. No point complicating the API IMO for something which is easy to by hand. |
|
Failing test seems to need |
|
Is this PR still waiting for the changes suggested in #4021 (comment)? By the way, I noticed that the code includes a workaround for a numpy 1.5.1 limitation. Now that support for numpy versions < 1.6 have been dropped, this workaround could be removed. |
|
Mostly. It still needs a little work on updating economic decomposition for Some code review would also be good. Even though there haven't been any On Wednesday, January 7, 2015, argriffing notifications@github.com wrote:
|
| for 1D arrays these are equivalent. This is Numpy's gh-2287, fixed in | ||
| 9b8ff38. | ||
|
|
||
| FIXME: Is it worth only applying this for numpy 1.5.1? |
There was a problem hiding this comment.
That commit is in numpy 1.6.2, which is now the lowest supported version. So this can go.
This is implemented by apply k consecutive rank-1 updates rather than attempting to block this algorithm. As part of this, the input validation code in qr_update was also reworked.
Works for me, but travis seems to fail depending on the exact configuration used.
This isn't so great yet for inserting a column vector that lies in the span of q.dot(r).
The output is not a QR decomposition if the column inserted lies in the span of Q. This will now fail in that case instead of returning an incorrect result.
7381603 to
da0616f
Compare
|
Rebased. |
|
Okay, 6 days since the message on the mailing list that this is about to be merged. Time to green-button it I think, then it gets to sit in master for a while before branching for the next release. |
ENH: add routines for updating QR decompositions
|
Or better, it could be used inside |
|
@rgommers Nice ! I'm not sure I'll have time to update place-poles with it before 0.16 though. |
|
No worries. If you may not be able to do that for a long time, it may be an idea to open an enhancement issue that explains what needs to be done. |
|
Should |
|
Indeed. Want to fix that, or should I? |
|
You can fix it while I let my scipy build recover from the numpy deprecation flag incompatibility :) |
|
done |
This pull request goes on top of #4021, accordingly only the last few commits are relevant here.
Here three new functions are added to perform various modifications of qr decompositions:
qr_update,qr_insertandqr_deletefor performing rank-k updates, adding row(s) or column(s) and removing row(s) or column(s).The methods used are essentially those in [1], with some input from [2].
I've endeavored to write this such that updates can be performed in-place as much as possible. This is, for instance, why the rank-1 update and the rank-k update algorithms are slightly different. This allows the rank-1 update to support updating
qandrwhen they are stored in C order. (In this spirit, it might be worth changinglinalg.qrto return fortran ordered 'r', since at present it always returns a C orderedr, and a F orderedq. Numpy's qr returns both as C order.)Edit: both are done (3/18/15)
Some outstanding things:
overwrite_qrshould probably default toFalseinstead of true.References