我在处理排名不足的情况时非常依赖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))
我们不能使用以下solve()反转此矩阵：
solve(badX)
## Error in solve.default(badX): system is computationally singular: reciprocal condition number = 5.55308e-18
然而，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
这是一个非常丑陋的结果。我尝试改变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")))
。我从中获得了其转置的 QR 分解，发现它应该是等级 4：
badQR <- qr(t(badL))
badQR$rank
## [1] 4
上述矩阵badX等于qr.R(badQR)[1:4, 1:4]，根据秩计算，它应该是一个全秩上三角矩阵。
我的补救措施似乎是使用zapsmall()以便我获得正确的排名......
qr(zapsmall(t(badL)))$rank
## [1] 1
我的问题是，为什么会这样？如果你看badL，很明显它有三个零行，只有第二行是非零的。我本以为qr()的枢轴方法会更好地处理这个。有没有更好的方法来获得更可靠的代码？
我正在运行 Windows 11 专业版，版本 10.0.22000 build 22000。这是我的 R 系统信息。
R.Version()
## $platform
## [1] "x86_64-w64-mingw32"
##
## $arch
## [1] "x86_64"
##
## $os
## [1] "mingw32"
##
## $crt
## [1] "ucrt"
##
## $system
## [1] "x86_64, mingw32"
##
## $status
## [1] ""
##
## $major
## [1] "4"
##
## $minor
## [1] "2.0"
##
## $year
## [1] "2022"
##
## $month
## [1] "04"
##
## $day
## [1] "22"
##
## $`svn rev`
## [1] "82229"
##
## $language
## [1] "R"
##
## $version.string
## [1] "R version 4.2.0 (2022-04-22 ucrt)"
##
## $nickname
## [1] "Vigorous Calisthenics"
创建于 2022-06-21 由 reprex 软件包 (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
。特别是(confounded)部分。此示例使用双因子模型，其中 2 个水平wool,3 个水平tension;但是，数据中省略了其中一个因子组合，这意味着对于tension主效应和wool:tension相互作用效应中的每一个，我们只能估计 1 d.f.，而对于wool则完全没有主效应。对于 5 个非空单元格的所有可能对比，有 4 d.f. 剩余的 2 d.f.，这些在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
。以及设计中所有单元格之间的对比：
> 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
我使用的方法是将tension和wool:tension的可估计函数结合起来，得到其转置的QR分解。然后我使用qr.resid()和上述单元格对比的转置。这给我们留下了(再次转置后)为(confounded)显示的可估计函数。该矩阵有 4 行，但其排名仅为 2，由该结果的 QR 分解确定;然后我提取 R 部分的 2x2 部分以完成F统计量的计算。
这个问题开头的例子类似，但模型更大、更复杂;badL矩阵是上述qr.resid()过程的结果。在这种情况下，其中一些行可以说应该是零。我目前的解决方法是检查R的对角线元素(OP中的badR)并选择超过绝对阈值的元素。
这里的基本思想是，我需要将所有对比的矩阵分解为两部分 - 已知的可估计函数和剩余函数。其中一个有趣的方面是，后半部分的等级是已知的，我没有利用这一事实。在未来的开发中，很可能是，根据@duffymo，使用SVD而不是这些带有qr.resid()的回旋。总是有新的东西需要学习...
你抱怨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
现在让我们看看你的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))
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
QR 分解在数值上是稳定的，因为对不同列的规模(或大小、幅度)不太敏感。 (其他因式分解，如 LU(solve基于它)和 SVD 可以。根据定义，这种分解确实
X = Q R
如果我们通过右乘全秩对角矩阵D来重新缩放X的列，则 QR 分解不会改变。
X D = Q R D
因此，让我们看一下应用 QR 分解的大矩阵t(badL)。我称之为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))
# [,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
让我们重新缩放其列，以便每列都有欧几里得范数(L2 范数，2 范数)1。
norm2 <- sqrt(colSums(X ^ 2))
XD <- X * rep(1 / norm2, each = nrow(X))
# [,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
你现在怎么看？它仍然是一个只有一个非零列的矩阵吗？尽管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
我们如何知道
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的条件数。
我建议你更喜欢奇异值分解。 如果矩阵排名不足，它将为您提供最佳解决方案。 下面是如何在 R 中使用它的示例。