lm()：LINPACK / LAPACK中QR分解返回的qraux是什么

Question

rich.main3

是 R 中的线性模型。我理解列表的其余元素，但我不明白

qraux

是什么。文档指出它是

长度为 ncol(x) 的向量，其中包含有关 old{Q}" 的附加信息。

这意味着什么附加信息？

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

想必你不知道QR分解是如何计算的。我用 LaTeX 编写了以下内容，可能会帮助您澄清这一点。当然，在编程网站上，我需要向您展示一些代码。最后我给大家提供一个R函数计算Householder反射的玩具。

---

住户反映矩阵

户主转型

户主 QR 分解（无旋转）

二维码的紧凑存储和缩放

---

LAPACK 辅助例程

dlarfg

正在执行 Householder 变换。我还编写了以下玩具 R 函数用于演示：

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
  }

假设我们有一个输入向量

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

我的函数给出：

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

Answer 2

希望以下内容能让一些学习者摆脱痛苦，因为我花了比我愿意承认的更多的时间来解决这个问题（还有一些令人讨厌的符号约定超出了我的兴趣）。所以就这样了：

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))
  rho = - sign(chi1) * norm
  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(t(augm)), as.vector(augm),'*')
  if(i < r){lwrT[(i + 1):r, i] <- H[(i + 1):r, i]}else{lwrT[i,i] <- tau}
  A <- H %*% A
}

qraux # As calculated above:
# [1] 1.8815031249 1.7072846053 1.9299080843 0.5442909928
qr(m)$qraux # Built-in function:
# [1] 1.8815031249 1.7072846053 1.9299080843 0.5442909928
qr(m)$qr # Built-in function:
              [,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
lwrT # As calculated above:
               [,1]          [,2]         [,3]         [,4]
 #[1,]  0.0000000000  0.0000000000 0.0000000000 0.0000000000
 #[2,] -0.4207603289  0.0000000000 0.0000000000 0.0000000000
 #[3,]  0.0969246757 -0.6907309238 0.0000000000 0.0000000000
 #[4,]  0.1910983876 -0.1504635437 0.3677919994 0.5442909928

来源这里：