对“qr()”的信心动摇

Question

对“qr()”的信心动摇

rvl*_*rvl 8 r matrix qr-decomposition matrix-decomposition

我在处理排名不足的情况时非常依赖该qr()函数，但最近遇到了一些它无法正常工作的示例。考虑下面的矩阵badX：

badX <-
structure(c(-1.641906809157e-10, 0, 0, 0, 0, -0.5, 0, 0, -1.10482935525559e-16, 
            0, -3.06266685765538e-17, 0, -4.83736007092039e-17, 0, -3.14414492582296e-18, 
            -3.06158275836099e-18), dim = c(4L, 4L), dimnames = list(c("(Intercept)", 
            "A2", "A3", "B2"), NULL))

Run Code Online (Sandbox Code Playgroud)

我们不能使用以下方法反转该矩阵solve()：

solve(badX)
## Error in solve.default(badX): system is computationally singular: reciprocal condition number = 5.55308e-18

Run Code Online (Sandbox Code Playgroud)

然而qr()，其相关例程认为该矩阵的秩为 4，并且可以对其求逆：

qr(badX)$rank
## [1] 4

qr.solve(badX)
##             [,1] [,2]          [,3]          [,4]
## [1,] -6090479645    0  2.197085e+10  7.366741e+10
## [2,]           0   -2  0.000000e+00  0.000000e+00
## [3,]           0    0 -3.265128e+16  3.353179e+16
## [4,]           0    0  0.000000e+00 -3.266284e+17

Run Code Online (Sandbox Code Playgroud)

这是一个相当难看的结果。我尝试改变tol论点，但结果没有改变。

对于上下文，此结果的来源是此对比矩阵：

badL <-
structure(c(0, 0, 0, 0, 0, -9.89189274870351e-11, 0, -5.55111512312578e-17, 
    -2.77555756156289e-17, 1.11022302462516e-16, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.25, 0, 0, 0, 0, -0.25, 0, 0, 
    0, 9.89189274870351e-11, 0, 5.55111512312578e-17, 2.77555756156289e-17, 
    -1.11022302462516e-16, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, -4.23939184015843e-11, 0, -4.16333634234434e-17, -1.38777878078145e-17, 
    5.55111512312578e-17, 0, 0, 0, 0, 0, -4.23939184015843e-11, 0, 
    -4.16333634234434e-17, -1.38777878078145e-17, 5.55111512312578e-17, 
    0, 0, 0, 0, 0, 0, 0.25, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.25, 0, 0, 
    0, 0, 0, 0, 0, 0, 4.23939184015843e-11, 0, 4.16333634234434e-17, 
    1.38777878078145e-17, -5.55111512312578e-17, 0, 0, 0, 0, 0, -1.41313127284714e-11, 
    0, -6.93889390390723e-18, -6.93889390390723e-18, 1.38777878078145e-17, 
    4.23939184015843e-11, 0, 4.16333634234434e-17, 1.38777878078145e-17, 
    -5.55111512312578e-17, 0, 0, 0, 0, 0), dim = c(5L, 24L), dimnames = list(
    NULL, c("(Intercept)", "A2", "A3", "B2", "B3", "C2", "C3", 
    "A2:B2", "A3:B2", "A2:B3", "A3:B3", "A2:C2", "A3:C2", "A2:C3", 
    "A3:C3", "B2:C2", "B3:C2", "B2:C3", "B3:C3", "A2:B2:C2", 
    "A3:B2:C2", "A3:B3:C2", "A2:B2:C3", "A3:B2:C3")))

Run Code Online (Sandbox Code Playgroud)

...从中我获得了其转置的 QR 分解，发现它应该是 4 级的：

badQR <- qr(t(badL))
badQR$rank
## [1] 4

Run Code Online (Sandbox Code Playgroud)

上述矩阵badX等于qr.R(badQR)[1:4, 1:4]基于秩计算的满秩上三角矩阵。

我的补救措施似乎是使用zapsmall()以便我获得正确的排名......

qr(zapsmall(t(badL)))$rank
## [1] 1

Run Code Online (Sandbox Code Playgroud)

我的问题是，为什么会发生这种情况？如果您查看badL，很明显它具有三个零行，并且只有第二行非零。我本以为qr()的旋转方法会更适合这个。有没有更好的方法来获得更可靠的代码？

我运行的是 Windows 11 Pro，版本 10.0.22000 build 22000。这是我的 R 系统信息。

solve(badX)
## Error in solve.default(badX): system is computationally singular: reciprocal condition number = 5.55308e-18

Run Code Online (Sandbox Code Playgroud)

^{由reprex 包于 2022 年 6 月 21 日创建(v2.0.1)}

有关上下文的更多信息

出现这个问题是因为我试图在emmeans包中产生这样的结果（对于一个更简单的例子）：

> (jt = joint_tests(warpx.emm))
 model term   df1 df2 F.ratio p.value note
 tension        1  37   5.741  0.0217    e
 wool:tension   1  37   5.867  0.0204    e
 (confounded)   2  37   7.008  0.0026  d e

d: df1 reduced due to linear dependence 
e: df1 reduced due to non-estimability

Run Code Online (Sandbox Code Playgroud)

...特别是这(confounded)部分。此示例使用具有wool2 个水平和tension3 个水平的双因素模型；然而，数据中省略了其中一个因子组合，这意味着我们只能估计tension主效应和wool:tension交互效应中每一个的 1 df，并且根本没有主效应wool。5 个非空单元格的所有可能对比有 4 个 df，还剩下 2 个 df，这些都在零件中confounded)。

计算基于相关的可估计函数：

> attr(jt, "est.fcns")
$tension
     (Intercept) woolB tensionM tensionH woolB:tensionM woolB:tensionH
[1,]           0     0        1        0            0.5              0

$`wool:tension`
     (Intercept) woolB tensionM tensionH woolB:tensionM woolB:tensionH
[1,]           0     0        0        0              1              0

$`(confounded)`
     (Intercept) woolB tensionM tensionH woolB:tensionM woolB:tensionH
[1,]           0    -1        0        0              0              0
[2,]           0     1        0        0              0              0
[3,]           0    -1        0        0              0              0
[4,]           0    -1        0        1              0              0

Run Code Online (Sandbox Code Playgroud)

...以及设计中所有单元之间的对比：

> contrast(warpx.emm, "consec")@linfct
     (Intercept) woolB tensionM tensionH woolB:tensionM woolB:tensionH
[1,]           0     1        0        0              0              0
[2,]           0    -1        1        0              0              0
[3,]           0     1        0        0              1              0
[4,]           0    -1       -1        1             -1              0
[5,]           0     1        0        0              0              1

Run Code Online (Sandbox Code Playgroud)

我使用的方法是结合tension和的可估计函数wool:tension，并获得其转置的QR分解。然后我将qr.resid()其与上述单元格的转置进行对比。这给我们留下了（再次转置之后）所示的可估计函数(confounded)。该矩阵有 4 行，但其秩仅为 2，这是由该结果的 QR 分解确定的；然后我提取R部分的2x2部分来完成F统计量的计算。

本问题开头的示例类似，但模型更大、更复杂；矩阵是上述过程badL的结果。qr.resid()在这种情况下，其中一些行可以说应该为零。我目前的解决方法是检查 R 的对角线元素（badR在 OP 中）并选择那些超过绝对阈值的元素。

这里的基本思想是，我需要将所有对比矩阵分解为两部分——已知的可估计函数和剩余部分。一个有趣的方面是后一部分的排名是已知的，但我没有利用这一事实。在未来的开发中，根据@duffymo，很可能使用 SVD 而不是这些带有qr.resid(). 总是有新的东西需要学习...

Answer 1

李哲源*_*李哲源 7

您抱怨solve无法反转似乎满秩的矩阵（根据qr）。你认为这solve是在做正确的事情，但事实qr并非如此。

好吧，不要相信solve。这不是一个稳健的数值过程，我们可以很容易地欺骗它。这是一个对角矩阵。它当然是可逆的（通过简单地反转其对角线元素），但solve就是做不到。

D <- diag(c(1, 1e-20))
#     [,1]  [,2]
#[1,]    1 0e+00
#[2,]    0 1e-20

solve(D)
#Error in solve.default(D) : 
#  system is computationally singular: reciprocal condition number = 1e-20

Dinv <- diag(c(1, 1e+20))

## an identity matrix, as expected
D %*% Dinv
#     [,1] [,2]
#[1,]    1    0
#[2,]    0    1

## an identity matrix, as expected
Dinv %*% D
#     [,1] [,2]
#[1,]    1    0
#[2,]    0    1

Run Code Online (Sandbox Code Playgroud)

现在让我们看看你的badX，我称之为R（因为它是 QR 分解返回的上三角矩阵）。

R <-
structure(c(-1.641906809157e-10, 0, 0, 0, 0, -0.5, 0, 0, -1.10482935525559e-16, 
            0, -3.06266685765538e-17, 0, -4.83736007092039e-17, 0, -3.14414492582296e-18, 
            -3.06158275836099e-18), dim = c(4L, 4L))

Run Code Online (Sandbox Code Playgroud)

solve不能反转它，但qr.solve给你一个适当的逆矩阵。

Rinv <- qr.solve(R)

## an identity matrix, as expected
R %*% Rinv
#     [,1] [,2] [,3]         [,4]
#[1,]    1    0    0 1.776357e-15
#[2,]    0    1    0 0.000000e+00
#[3,]    0    0    1 0.000000e+00
#[4,]    0    0    0 1.000000e+00

## an identity matrix, as expected
Rinv %*% R
#     [,1] [,2] [,3]         [,4]
#[1,]    1    0    0 5.293956e-23
#[2,]    0    1    0 0.000000e+00
#[3,]    0    0    1 1.387779e-17
#[4,]    0    0    0 1.000000e+00

Run Code Online (Sandbox Code Playgroud)

QR 分解在数值上是稳定的，因为它对不同列的尺度（或大小、大小）不太敏感。（其他分解，如 LU（基于其solve）和 SVD 也是如此。）根据定义，此分解不会

X=QR

如果我们通过右乘满秩对角矩阵D来重新缩放X的列，则 QR 分解不会改变。

XD = QRD

让我们看看t(badL)应用 QR 分解的大矩阵。我称之为X。

X <- structure(c(0, -9.89189274870351e-11, 0, 0, 0, 0, 0, 9.89189274870351e-11, 0, 0, 0, -4.23939184015843e-11, 0, -4.23939184015843e-11, 0, 0, 0, 0, 0, 4.23939184015843e-11, 0, -1.41313127284714e-11, 4.23939184015843e-11, 0, 0, 0, 0, 0, 0, -0.25, -0.25, 0, 0, 0, 0, 0, 0, 0, 0, 0.25, 0, 0.25, 0, 0, 0, 0, 0, 0, 0, -5.55111512312578e-17, 0, 0, 0, 0, 0, 5.55111512312578e-17, 0, 0, 0, -4.16333634234434e-17, 0, -4.16333634234434e-17, 0, 0, 0, 0, 0, 4.16333634234434e-17, 0, -6.93889390390723e-18, 4.16333634234434e-17, 0, 0, -2.77555756156289e-17, 0, 0, 0, 0, 0, 2.77555756156289e-17, 0, 0, 0, -1.38777878078145e-17, 0, -1.38777878078145e-17, 0, 0, 0, 0, 0, 1.38777878078145e-17, 0, -6.93889390390723e-18, 1.38777878078145e-17, 0, 0, 1.11022302462516e-16, 0, 0, 0, 0, 0, -1.11022302462516e-16, 0, 0, 0, 5.55111512312578e-17, 0, 5.55111512312578e-17, 0, 0, 0, 0, 0, -5.55111512312578e-17, 0, 1.38777878078145e-17, -5.55111512312578e-17, 0), dim = c(24L, 5L))
Run Code Online (Sandbox Code Playgroud)
# [,1] [,2] [,3] [,4] [,5] # [1,] 0.000000e+00 0.00 0.000000e+00 0.000000e+00 0.000000e+00 # [2,] -9.891893e-11 0.00 -5.551115e-17 -2.775558e-17 1.110223e-16 # [3,] 0.000000e+00 0.00 0.000000e+00 0.000000e+00 0.000000e+00 # [4,] 0.000000e+00 0.00 0.000000e+00 0.000000e+00 0.000000e+00 # [5,] 0.000000e+00 0.00 0.000000e+00 0.000000e+00 0.000000e+00 # [6,] 0.000000e+00 -0.25 0.000000e+00 0.000000e+00 0.000000e+00 # [7,] 0.000000e+00 -0.25 0.000000e+00 0.000000e+00 0.000000e+00 # [8,] 9.891893e-11 0.00 5.551115e-17 2.775558e-17 -1.110223e-16 # [9,] 0.000000e+00 0.00 0.000000e+00 0.000000e+00 0.000000e+00 #[10,] 0.000000e+00 0.00 0.000000e+00 0.000000e+00 0.000000e+00 #[11,] 0.000000e+00 0.00 0.000000e+00 0.000000e+00 0.000000e+00 #[12,] -4.239392e-11 0.00 -4.163336e-17 -1.387779e-17 5.551115e-17 #[13,] 0.000000e+00 0.00 0.000000e+00 0.000000e+00 0.000000e+00 #[14,] -4.239392e-11 0.00 -4.163336e-17 -1.387779e-17 5.551115e-17 #[15,] 0.000000e+00 0.00 0.000000e+00 0.000000e+00 0.000000e+00 #[16,] 0.000000e+00 0.25 0.000000e+00 0.000000e+00 0.000000e+00 #[17,] 0.000000e+00 0.00 0.000000e+00 0.000000e+00 0.000000e+00 #[18,] 0.000000e+00 0.25 0.000000e+00 0.000000e+00 0.000000e+00 #[19,] 0.000000e+00 0.00 0.000000e+00 0.000000e+00 0.000000e+00 #[20,] 4.239392e-11 0.00 4.163336e-17 1.387779e-17 -5.551115e-17 #[21,] 0.000000e+00 0.00 0.000000e+00 0.000000e+00 0.000000e+00 #[22,] -1.413131e-11 0.00 -6.938894e-18 -6.938894e-18 1.387779e-17 #[23,] 4.239392e-11 0.00 4.163336e-17 1.387779e-17 -5.551115e-17 #[24,] 0.000000e+00 0.00 0.000000e+00 0.000000e+00 0.000000e+00
Run Code Online (Sandbox Code Playgroud)
让我们重新缩放其列，以便每列都具有欧几里德范数（L2 范数，2-范数）1。

norm2 <- sqrt(colSums(X ^ 2)) XD <- X * rep(1 / norm2, each = nrow(X))
Run Code Online (Sandbox Code Playgroud)
# [,1] [,2] [,3] [,4] [,5] # [1,] 0.00000000 0.0 0.00000000 0.0000000 0.00000000 # [2,] -0.60246371 0.0 -0.48418203 -0.5714286 0.57585260 # [3,] 0.00000000 0.0 0.00000000 0.0000000 0.00000000 # [4,] 0.00000000 0.0 0.00000000 0.0000000 0.00000000 # [5,] 0.00000000 0.0 0.00000000 0.0000000 0.00000000 # [6,] 0.00000000 -0.5 0.00000000 0.0000000 0.00000000 # [7,] 0.00000000 -0.5 0.00000000 0.0000000 0.00000000 # [8,] 0.60246371 0.0 0.48418203 0.5714286 -0.57585260 # [9,] 0.00000000 0.0 0.00000000 0.0000000 0.00000000 #[10,] 0.00000000 0.0 0.00000000 0.0000000 0.00000000 #[11,] 0.00000000 0.0 0.00000000 0.0000000 0.00000000 #[12,] -0.25819930 0.0 -0.36313652 -0.2857143 0.28792630 #[13,] 0.00000000 0.0 0.00000000 0.0000000 0.00000000 #[14,] -0.25819930 0.0 -0.36313652 -0.2857143 0.28792630 #[15,] 0.00000000 0.0 0.00000000 0.0000000 0.00000000 #[16,] 0.00000000 0.5 0.00000000 0.0000000 0.00000000 #[17,] 0.00000000 0.0 0.00000000 0.0000000 0.00000000 #[18,] 0.00000000 0.5 0.00000000 0.0000000 0.00000000 #[19,] 0.00000000 0.0 0.00000000 0.0000000 0.00000000 #[20,] 0.25819930 0.0 0.36313652 0.2857143 -0.28792630 #[21,] 0.00000000 0.0 0.00000000 0.0000000 0.00000000 #[22,] -0.08606647 0.0 -0.06052275 -0.1428571 0.07198158 #[23,] 0.25819930 0.0 0.36313652 0.2857143 -0.28792630 #[24,] 0.00000000 0.0 0.00000000 0.0000000 0.00000000
Run Code Online (Sandbox Code Playgroud)
现在你觉得怎么样？它仍然是一个只有一个非零列的矩阵吗？虽然qr(X)在 QR 分解之前实际上并不首先重新缩放所有列，但查看XD确实可以帮助您更好地理解为什么 QR 分解更稳健。

如果您确实想干预，请不要使用zapsmall; 相反，阈值列的 2-范数。

X0 <- X X0[, norm2 < max(norm2) * sqrt(.Machine$double.eps)] <- 0 QR0 <- qr(X0) QR0$rank # [1] 1
Run Code Online (Sandbox Code Playgroud)

我们如何知道这sqrt(.Machine$double.eps)是一个合适的阈值？

sqrt(.Machine$double.eps)（大约 1e-8）和之间的任何阈值.Machine$double.eps (about 1e-16)都是合理的。使用.Machine$double.eps可恢复常规 QR 结果，为您提供排名 4。

“sqrt”阈值来自我们想要查看的情况X'X，它是的条件数的平方X。

归档时间：	3 年，2 月前
查看次数：	189 次
最近记录：	3 年，2 月前