2

我有一种情况，我有一个秩不足的 X 矩阵，但仍然需要回归系数。具体来说，我想要由提供的表现良好的系数base::qr

require(Matrix)
y <- c(230, 192, 195, 180, 200, 185, 0)
X <- new("dgCMatrix", i = c(0L, 1L, 2L, 3L, 4L, 5L, 6L, 0L, 1L, 2L, 
         3L, 4L, 5L, 1L, 3L, 5L, 2L, 3L, 4L, 5L, 1L, 0L), p = c(0L, 7L, 
         13L, 16L, 20L, 20L, 21L, 22L), Dim = c(7L, 7L), Dimnames = list(
         c("419", "420", "421", "422", "423", "424", ""), c("", "(Intercept)", 
         "A", "B1", "B2", "B3", "B4")), x = c(0.451764131343143, 
         0.434020584428936, 0.451764131343143, 0.434020584428936, 0.451764131343143, 
         0.451764131343143, 0.788669766125214, 0.788669766125214, 0.757693869580585, 
         0.788669766125214, 0.757693869580585, 0.788669766125214, 0.788669766125214, 
         0.757693869580585, 0.757693869580585, 0.788669766125214, 0.788669766125214, 
         0.757693869580585, 0.788669766125214, 0.788669766125214, 0.757693869580585, 
         0.788669766125214), factors = list())

qrB <- base::qr(X)
qrM <- Matrix::qr(X)
qr.coef(qrB, y)
#                 (Intercept)             A            B1            B2            B3            B4 
# -8.438370e-14  2.916303e+02 -1.441398e+01 -4.120863e+01            NA -2.381583e+01            NA 
qr.coef(qrM, y)
# [1] -9.730332e-01 -1.314143e+16 -1.597833e+01  1.314143e+16  0.000000e+00  1.314143e+16  1.314143e+16
# Warning message:
# In lengths(res@Dimnames) :
#   sparseQR_coef(): structurally rank deficient case: possibly WRONG zeros

因此，双方都认识到他们无法估计 B2 和估计 na NA ( base) 或零 ( Matrix)，我对任何一个结果都很好。但是，base::qr也认识到它需要删除第二列并选择最后一列（目前尚不清楚为什么，但它有效），同时Matrix::qr大胆地进行并估计一个类似于.Machine$double.eps.

矩阵也不太稳定（请原谅​​我，我忘记了哪种），因为平方和比基数更小

coef <- qr.coef(qrB, y)
coef[is.na(coef)] <- 0
sum((X %*% coef  - y)^2)
# [1] 15.17039
sum((X %*% qr.coef(qrM, y)  - y)^2)
# [1] 18.58891

有没有办法获得一种超级快速的方法，就像Matrix::qr工作一样base::qr并给出快速但稳定的结果？我意识到这并不完全是一个“大矩阵”，部分原因是它是一个 MWE。

Ben Bolker's answer here和两个 answers here谈论 Lapack 和 Linpack 之间的区别，但这些方法在base::qr. 我正在寻找一种方法来从Matrix::qr我可以使用的系数中获得修改后的结果。

4

0 回答 0