I have a matrix $A \in \mathbb{R}^{m \times n}$ where $m \gg n$ and I want to compute the full QR decomposition $A = QR$. Where $Q$ is an orthogonal $m \times m$ matrix.

Bishof & Van Loan (1987) describe that $Q$ can be represented as the sum of a low-rank matrix plus the identity matrix:

$$Q = I + W Y$$

Where $W \in \mathbb{R}^{m \times n}$ and $Y \in \mathbb{R}^{n \times m}$. I implemented this myself in Python, and it works! But it is still pretty slow because it involves a for loop over the $n$ columns. I'm curious if anyone has pointers for libraries that already implement this... If not I might try re-implementing in Julia.


A later paper, Schrieber and Van Loan 1989, describe a "compact" WY representation which is more storage efficient, $$ Q = I - Y T Y^T $$ where $T$ is upper triangular. The paper by Elmroth and Gustavson 2000 then use this compact WY representation to derive a recursive QR algorithm, which is quite elegant and efficient as it is based on Level-3 BLAS operations.

A C implementation of the Elmroth and Gustavson 2000 algorithm is in the GSL library: here.

LAPACK provides a fortran implemenation of Elmroth and Gustavson in the DGEQRT3 routine. The standard LAPACK QR routine (DGEQRF) also uses the WY approach, as does the rank revealing QR algorithm (technical note here). However these (non-recursive) algorithms are quite a bit more complex, and not so elegant.

Schreiber, Robert, and Charles Van Loan. "A storage-efficient WY representation for products of Householder transformations." SIAM Journal on Scientific and Statistical Computing 10.1 (1989): 53-57.

Elmroth, E. and Gustavson, F.G., 2000. Applying recursion to serial and parallel QR factorization leads to better performance. IBM Journal of Research and Development, 44(4), pp.605-624.


Just building on @vibe's answer. Here's a quick sketch of how to access and work with LAPACK to do this with Python/scipy. This uses the Elmroth & Gustavson algorithm under the hood.

from scipy.linalg.lapack import dgeqrt

def tall_qr(A):
    m, n = A.shape
    assert m > n
    V, T, INFO = dgeqrt(n, A)

    # Extract R from V.
    R = np.zeros_like(V)
    R[np.triu_indices_from(R)] = V[np.triu_indices_from(V)]

    # Extract elementary reflectors in V.
    V[np.arange(n), np.arange(n)] = 1.0
    V[np.triu_indices_from(V, 1)] = 0.0

    # To recover the dense orthogonal matrix use:
    #   Q = np.eye(m) - V @ T @ V.T
    # However, the dense matrix multiply (Q @ B) should be
    # slower than B - (V @ T) @ (V.T @ B) for large matrices.

    return (V, T), R
  • 1
    $\begingroup$ If you need to form the full $Q$ (which is rare and probably unnecessary), there is a more efficient way to do it (see the gsl_linalg_QR_unpack_r routine in the GSL link above). You correctly note that computing $Q B$ or $Q^T B$ is more efficiently done by applying the matrices $V,T,V^T$ individually to $B$. See the QTmat routine in the GSL link. Finally, note that the DGEQRF routine handles the case $N < M$ while E&G does not. It may be possible to extend E&G to this case, but I have not tried (and I'm not aware of others doing it either). $\endgroup$ – vibe Oct 24 '20 at 1:50
  • 1
    $\begingroup$ Also, if your matrix happens to be so tall that you cannot fit the entire matrix in memory at one time, you can use the TSQR algorithm (arxiv.org/abs/0806.2159), which can use E&G for the individual QR updates of a block of rows. $\endgroup$ – vibe Oct 24 '20 at 2:07
  • $\begingroup$ @vibe -- I'm trying to figure out if there's a blas / lapack routine that will do the low rank update B - (V @ T) @ (V.T @ B) inplace on B? Right now I'm just writing it out in numpy, but B has quite a large number of columns, so the sub-expression (V @ T) @ (V.T @ B) instantiates a large matrix in memory. DLARFB seems to be the right method to use? Unfortunately, it doesn't seem exposed in scipy.linalg.lapack... $\endgroup$ – ahwillia Oct 24 '20 at 4:40
  • $\begingroup$ You might try DGEMQRT, which looks like a wrapper function for DLARFB. There is also a DGEMQR function which looks like a wrapper for DGEMQRT. Also, is there a reason you don't want to use the standard DGEQRF function for your $QR$ decomposition? $\endgroup$ – vibe Oct 24 '20 at 6:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.