Julia: Rank 1 updates to Cholesky and QR factorizations

[ViralBShah]

Matlab has cholupdate for rank 1 updates and downdates to Cholesky factorizations. It would be nice to have something similar, but perhaps with better syntax.

LINPACK had these, but LAPACK left out the Cholesky update/downdate capabilities.

http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=2646

For the dense case, perhaps these should be implemented in julia. For the sparse case, I believe that CHOLMOD does provide these capabilities.

Updated by @andreasnoack:

• [x] Dense Cholesky
• [x] Sparse Cholesky
• [ ] Dense QR (actually there are three different storage schemes)

[ViralBShah]

Cc: @dmbates @andreasnoack

[mlubin]

I think qrupdate is the best package out there for this. It's GPL. Implementing this in Julia won't be a walk in the park; it's over 6000 lines of fortran (some of which is duplicated logic for different types of matrices).

[GordStephen]

Did this ever make it in? I can't find it anywhere, but I also don't know my way around linalg that well.

[ViralBShah]

This is still open.

[GordStephen]

Ok, good to know. I'm not really qualified to be adding cross-platform binary dependencies to Julia, but I did get a rudimentary version of cholupdate! (calling qrupdate) up and running for personal use on my local machine, if it's of use to anyone in the future: https://gist.github.com/GordStephen/be274faa22417d6f26a1

[ViralBShah]

Perhaps make it a package?

[GordStephen]

Good idea - I might not get to it this week but when I do I'll add in functionality for the other input types supported by qrupdate as well.

[andreasnoack]

We should be able to add pure Julia Cholesky updates without too much effort. We already have the necessary components and the algorithm is fairly simple. I can take a look. The QR row updates are probably more complicated.

[tkelman]

See @mpf's https://github.com/mpf/QRupdate.jl

[klkeys]

Any progress on sparse rank-1 Cholesky updates? I ask because I noticed that CHOLMOD is supposed to export them. See the reference publication here.

[tkelman]

If those are available in cholmod, I don't think we have wrappers for them yet.

[klkeys]

ahahaha sorry @tkelman I was updating the comment a minute too late. Are the wrappers difficult to code? I suddenly have occasion to use them, but I am not very good with C/Fortran...

[tkelman]

I get an access denied on that link. Writing ccall wrappers is conceptually not too bad, though CHOLMOD has a pretty confusing API at times. Some quality time with the API docs and/or headers along with http://docs.julialang.org/en/release-0.4/manual/calling-c-and-fortran-code/ and looking at the existing bindings in https://github.com/JuliaLang/julia/blob/master/base/sparse/cholmod.jl should get you started.

[papamarkou]

I need to use rank 1 Cholesky updates for an adaptive algorithm I am working on, to bring down complexity from O(n^3) to O(n^2), and was a bit unclear from this open issue if there is some preliminary yet operational support (either built-in or by invoking QRupdate)?

[bicycle1885]

LinAlg.lowrankupdate! and LinAlg.lowrankdowndate! will do. These are included in Julia 0.5.

[papamarkou]

Thanks, @bicycle1885.

[mfalt]

I implemented a wrapper to Cholmod for low rank updates of Sparse Cholesky in https://github.com/JuliaLang/julia/pull/21426 that solves the second checkbox.
Should we also consider exporting lowrankupdate(!) and lowrankdowndate(!)?

[andreasnoack]

Should we also consider exporting lowrankupdate(!) and lowrankdowndate(!)?

Yes. We should either export these from Base or more consistently not export LinAlg functions and instead require using LinAlg. For now, I think it would be fine to export them.

[ViralBShah]

@andreasnoack I thought I saw some related work for this but can't quite figure out where it happened. Can you chime in since you may know the latest on this?

[andreasnoack]

The sparse Cholesky was added in #21426. I've updated the top post. Was that what you were thinking of?

[pkofod]

Are QR updates hard to implement (wrap?), or is it just a rare feature to request? I'm asking because I'd love this functionality, but have no idea where the work would start or end.

[Nosferican]

I would also love being able to update a QR factorization by adding or deleting columns. Deleting columns is trivial by removing the columns from R, but missing being able to add columns. Don't know if applicable, but if an in-place version can be developed for replacing columns that would be great too.

This reference might help: Updating the QR factorization and the least squares problem

[andreasnoack]

Updating a QR factorization is not that hard in theory but, in contrast to the Cholesky factorization, the storage scheme required doesn't align with the way we store the QR factorizations. I can see from the documentation of qrupdate that Q is stored as a dense matrix in contrast to storing Q in Householder form which we do. In addition, if you add a column the sizes Q and R changes which is not permitted in our current version of QR. It might be possible to come up with a new factorization type, say QRUpdate, that allocated sufficient memory for the updates. That work would be a good candidate for a package though.

[Nosferican]

Would the same issues apply to replacing or dropping columns? I believe dropping columns could be easily implemented by zeroing out the column in R after computing it from the Householder operations. Replacing columns might be even more useful than adding. In my probably poor understanding, QR decomposition is the way to finding linearly independent / dependent columns for rectangular matrices (singular values and pivoted Cholesky don't give the location of the dependent columns; rref is the alternative, but equally expensive). In terms of least squares problems one may have all possible variables and then take a linear independent subset of these (dropping collinear columns). If one uses instrumental variable estimations, then one needs to replace the endogenous variables with their instrumented versions (replacement).

[mpf]

In case it helps, the package QRupdate.jl does Q-less QR updates: add/delete columns and add rows.

[Nosferican]

Thanks, I use both Q and R which makes deleting columns trivial, but adding columns was a feature I was missing. Absent a replacement procedure this seems like the best available solution.

[Nosferican]

After some preliminary benchmarking it seems the best combination is to use qrfact(hcat(X, Z)) and then update R using qrdelcol when there are a few linearly dependent columns... A more developed benchmark could identify when is better to refactor using qrfact based on m, n, and how many columns are dependent. Likewise, to add columns...

If I am not mistaken, it seems updating R is number of column updates complex while recomputing new R is number of rows and number of final columns complex... Since in least squares m tends to overly dominate n for most cases it would be best to update.

@mpf do you have some benchmarks or guidelines to advise when to use which method?
Also, would you recommend qrfact(::AbstractMatrix)[:R]or would it be possible to make a more efficient version in the package for full Q-less QR decomposition rather than building it column by column? I believe qrfact uses the projection multiplications for calculating R.

[pkofod]

Would you mind posting the script?

[Nosferican]

This is a very preliminary skeleton for a more comprehensive benchmark gist. Three functions are tested (1) building R using QRupdate, (2) updating R using QDupdate, (3) updating R by factoring R again. Each returns the linearly independent subset of two matrices, the gram matrix, and an indicator vector of which columns were independent.

If y'all develop a comprehensive benchmark ping me to see the results. I am interested in routines in least squares estimation.

[ViralBShah]

Should we rename this issue for the dense case? Or do we already have the dense QR update?

[pkofod]

Or do we already have the dense QR update?

We don't have dense QR updates I believe (@andreasnoack and I discussed this a few weeks ago in the context of quadratic programming, so something should have happened since then).

[Nosferican]

I believe there isn't any dense QR updates. As part of the implementation, I proposed having an option for using the updates if possible to compute the QR decomposition of the full-rank version of a matrix directly.