lm(): What is qraux returned by QR decomposition in LINPACK / LAPACK

Question

rich.main3 is a linear model in R. I understand the rest of the elements of the list but I don't get what qraux is. The documentation states that it is

a vector of length ncol(x) which contains additional information on \bold{Q}".

What additional information does it mean?

str(rich.main3$qr)

qr   : num [1:164, 1:147] -12.8062 0.0781 0.0781 0.0781 0.0781 ...


..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:164] "1" "2" "3" "4" ...
  .. ..$ : chr [1:147] "(Intercept)" "S2" "S3" "x1" ...
  ..- attr(*, "assign")= int [1:147] 0 1 1 2 3 4 5 6 7 8 ...
  ..- attr(*, "contrasts")=List of 3
  .. ..$ S    : chr "contr.treatment"
  .. ..$ ID   : chr "contr.treatment"
  .. ..$ Block: chr "contr.treatment"
 $ qraux: num [1:147] 1.08 1.06 1.16 1.21 1.27 ...
 $ pivot: int [1:147] 1 2 3 4 5 6 7 8 10 11 ...
 $ tol  : num 1e-07
 $ rank : int 21
 - attr(*, "class")= chr "qr"

Answer 1

Hopefully the following can put some learners out of their misery since it took me more hours than I care to admit to sort this out. So here it goes...

The mystery behind the vector with "additional information" and the "compact form" storage hopefully clarified.

qr a matrix with the same dimensions as x. The upper triangle contains the R of the decomposition and the lower triangle contains information on the Q of the decomposition (stored in compact form). Note that the storage used by DQRDC (LINPACK) and DGEQP3 (LAPACK) differs.

qraux a vector of length ncol(x) which contains additional information on Q.

Source in here.

Here is an illustration of the different "leading" column vectors (they will change with each update of A), and the equations as presented in the quoted source above for tau (and the qraux vector):

set.seed(2024)
options(digits=10)

# Since it's about understanding the convention of the qr()
# I'm sticking to square matrices. s is the number of rows and cols:
s = 4
m = matrix(rnorm(s^2), nrow = s)

# chi1 is the leading entry of each col vector in the matrix m (aka as A)
# keeping in mind that the matrix will be updated in a loop
# by multiplying it on the left by the Householder matrices:
chi1 = m[1,1] # This is the initial entry out of the original matrix.
# The initial col vector is precisely the first col in the orig. mat.
# And its 2-norm is:
norm = sqrt(sum(m[,1]^2))
# The Householder vector will be constructed by taking the entries of
# the column matrix at each iteration of the products H H H A
# below the leading entry in the decreasing sub-matrices (chi1):
x2 = m[2:nrow(m), 1] # This is the initial one.
# nu1 is the difference between the original vector and the mirror image
# avoiding catastrophic cancellation, ie chi1 + norm(x) * e_1:
nu1 = chi1 + sign(chi1) * norm
# nu1 will standardize u2:
u2 = x2 / nu1
# Let's see this applied to this first column vector in the matrix A
(tau = 2 / (1 + t(u2)%*%u2))

# Now let's get all the values of the output in qraux and in the lower
# triangular matrix of the output:
# Initializing a vector to collect values below the diag
# They will be output with a sign difference wrt to the R base function.
lwrT = matrix(0, s, s)
# Initializing an empty vector to collect the values of qraux:
qraux <- numeric()
A = m
r = nrow(m)
for (i in 1:r){
  chi1 = A[i,i]
  norm = sqrt(sum(A[i:r, i]^2))
  if(i < r){x2 = as.matrix(A[(i + 1):r, i])}else{x2 = 0}
  nu1 = chi1 + sign(chi1) * norm
  u2 = x2 / nu1
  if(i < r){tau = 2 / (1 + t(u2) %*% u2)}else{tau = sign(chi1) * nu1 / 2}
  qraux[i] = tau
  if(i < r) {augm = c(rep(0, i - 1), 1, u2)}
          else
            {augm = c(rep(0, i - 1), 0)}
  H <- diag(r) - as.numeric(tau) * 
    outer(as.vector(augm), as.vector(t(augm)),'*')
  if(i < r){lwrT[(i + 1):r, i] <- - H[(i + 1):r, i]}else{}
  A <- H %*% A
}
# After all iterations multiplying H A we end up producing
# R or the upper triangular (including the diagonal in qr()):
round(A,9)
 #       [,1]         [,2]        [,3]         [,4]
 #[1,] -1.11397156 -1.537094686 1.479110318 -1.202532687
 #[2,]  0.00000000 -0.975037638 1.587258014  2.295022409
 #[3,]  0.00000000  0.000000000 1.032380986  2.124104691
 #[4,]  0.00000000  0.000000000 0.000000000 -0.544290993
# The lower triangular containing essentially tau * uTu:
lwrT
 #         [,1]         [,2]          [,3] [,4]
 #[1,]  0.0000000000 0.0000000000  0.0000000000    0
 #[2,]  0.4207603289 0.0000000000  0.0000000000    0
 #[3,] -0.0969246757 0.6907309238  0.0000000000    0
 #[4,] -0.1910983876 0.1504635437 -0.3677919994    0
# They fit perfectly to replicate the output of qr():
qr(m)$qr
#          [,1]          [,2]          [,3]          [,4]
#[1,] -1.1139715602 -1.5370946856  1.4791103176 -1.2025326872
#[2,]  0.4207603289 -0.9750376385  1.5872580139  2.2950224086
#[3,] -0.0969246757  0.6907309238  1.0323809861  2.1241046912
#[4,] -0.1910983876  0.1504635437 -0.3677919994 -0.5442909928
# Likewise the output of the vector in the calculations above
# for tau:
qraux
#[1] 1.8815031249 1.7072846053 1.9299080843 0.5442909928
# replicate the built-in function qraux:
qr(m)$qraux
#[1] 1.8815031249 1.7072846053 1.9299080843 0.5442909928

Answer 2

Presumably you don't know how QR factorization is computed. I wrote the following in LaTeX which might help you clarify this. Surely on a programming site I need to show you some code. In the end I offer you a toy R function computing Householder reflection.

---

Householder reflection matrix

Householder transformation

Householder QR factorization (without pivoting)

Compact storage of QR and rescaling

---

The LAPACK auxiliary routine dlarfg is performing Householder transform. I have also written the following toy R function for demonstration:

dlarfg <- function (x) {
  beta <- -1 * sign(x[1]) * sqrt(as.numeric(crossprod(x)))
  v <- c(1, x[-1] / (x[1] - beta))
  tau <- 1 - x[1] / beta
  y <- c(beta, rep(0, length(x)-1L))
  packed_yv <- c(beta, v[-1])
  oo <- cbind(x, y, v, packed_yv)
  attr(oo, "tau") <- tau
  oo
  }

Suppose we have an input vector

set.seed(0); x <- rnorm(5)

my function gives:

dlarfg(x)
#              x         y           v   packed_yv
#[1,]  1.2629543 -2.293655  1.00000000 -2.29365466
#[2,] -0.3262334  0.000000 -0.09172596 -0.09172596
#[3,]  1.3297993  0.000000  0.37389527  0.37389527
#[4,]  1.2724293  0.000000  0.35776475  0.35776475
#[5,]  0.4146414  0.000000  0.11658336  0.11658336
#attr(,"tau")
#[1] 1.55063