How to most efficiently use QR-decomposition in Julia?

Question

Avoiding array allocations is good for performance. However, I have yet to understand what is the most possible efficient way one can perform a QR decomposition of a matrix A. (note: both Q and R matrices are needed)

Simply using Q, R = qr(A) is probably not the best idea, since it allocates both Q and R, where both could be re-allocated.

The function qrfact allows one to store factorization in a packed format. However, I would still write afterwards: F = qrfact(A); Q = F[:Q]; R = F[:R] once again allocating new arrays for Q and R. Finally, the documentation also suggests the qrfact! function, which saves space by overwriting the input A, instead of creating a copy. However, if one uses F = qrfact!(A) the over-written A is not useful in the sense that it is not either Q or R, which one (specifically, I) would need.

So my two questions are:

1. What is the best/most efficient way to perform a QR decomposition if you only care about the matrices Q and R and you have no problem re-allocating them.

2. What is actually written in the matrix A when one calls qrfact!(A) ?

Cloninger

Why do you need to do Q = F[:Q]; R = F[:R] afterwards? Can't you just use F[:Q] as it is? – Enmesh

There are references to the storage method in the Julia documentation. Specifically this paper was linked to. But in fact, you would want to look at this julia code – Blackout

@MichaelK.Borregaard It doesn't matter for me if I use Q or F[:Q]. The issue is how to avoid the allocation done from F[:Q] (if any is done). – Cloninger

Oh damn, unfortunately I need to operate on the F[:Q] matrix because I need to change sign on some of the columns (qr-decomposition to find lyapunov exponents) – Cloninger

I have some serious doubts that the cost to allocate an array (effectively O(1)) is even detectable in an algorithm that involves performing the O(n^3) QR decomposition. Have you profiled your code and verified that you actually need to avoid these allocations? – Moa

Answer 1

In

F = qrfact!(A)

or

F = qrfact(A)

F[:Q] and F[:R] do not allocate new dense arrays; they are simply views over the packed format from which Q and R are easily computed. This means that qrfact!(A) doesn't need to allocate arrays for Q and R, it simply computes the packed format in place for A.

However, that also means that F[:Q] and F[:R] cannot be mutated. If you need to modify one of them for whatever reason, you will need to collect it into a mutable Array, and this will certainly allocate. It will still be more efficient to use qrfact!(A) instead of qrfact(A), because the latter will allocate space for the packed QR factorization as well as for the collected Array.

Moa

Thanks for the answer. After looking around more I also came to the same conclusion (F[:Q] doesn't allocate). I've also understood what qrfact!(A) does to A: It has the same elements as R, but it is not upper-triagonal, so I don't know whats going on with the lower triangle though. P.s.: No I have not found this to be the reason my code is slow. It was simply an opportunity for me to get a better understanding of how to use qr with Julia. – Cloninger

The format we use to store the QR factorization is described in netlib.org/lapack/explore-html/dd/d9a/… so basically, Q is stored in the lower triangle and the extra T matrix. – Drin