7

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

George Datseris
  • 411
  • 4
  • 14
  • Why do you need to do `Q = F[:Q]; R = F[:R]` afterwards? Can't you just use `F[:Q]` as it is? – Michael K. Borregaard May 28 '17 at 12:04
  • There are references to the storage method in the Julia documentation. Specifically this [paper](https://www.cs.cornell.edu/cv/ResearchPDF/A%20Storage-Efficient%20WY%20Representation%20for%20Products%20of%20Householder%20Transformations.pdf) was linked to. But in fact, you would want to look at this [julia code](https://github.com/JuliaLang/julia/blob/82e96382251a1e163860df03c846acec86d07c36/base/linalg/qr.jl#L43-L87) – Dan Getz May 28 '17 at 12:52
  • @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). – George Datseris May 28 '17 at 13:13
  • 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) – George Datseris May 28 '17 at 14:20
  • 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? – Fengyang Wang May 28 '17 at 17:07

1 Answers1

6

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.

Fengyang Wang
  • 11,901
  • 2
  • 38
  • 67
  • 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. – George Datseris May 29 '17 at 13:22
  • The format we use to store the QR factorization is described in http://www.netlib.org/lapack/explore-html/dd/d9a/group__double_g_ecomputational_ga3ad112f2b0890b3815e696628906f30c.html#ga3ad112f2b0890b3815e696628906f30c so basically, `Q` is stored in the lower triangle and the extra `T` matrix. – Andreas Noack May 30 '17 at 01:51