matrix r qr-decomposition matrix-decomposition

动摇了对'qr（）'的信仰

发布于 2025-02-09 07:06:45 字数 5988 浏览 2 评论 0原文

我在处理排名不足的情况方面非常依赖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 Pro，版本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"

^由

创建在上下文上，

这个问题出现了，因为我试图在 ememeans 软件包中产生这样的结果（以简单的示例）：

> (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

...尤其是（混淆）部分。此示例是带有羊毛在2级的两因素模型，张力在3级;但是，数据中省略了因子组合之一，这意味着我们只能估算每个张力主要效果和羊毛的1 df：张力交互效果，对于羊毛，根本没有主要效果。对于5个非空单元的所有可能对比度，有4个df，剩下2个df，而这些df在混淆中）部分。

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

> 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

我使用的方法是结合张力>张力和羊毛的可估计功能，并获得其转置的QR分解。然后，我将其使用QR.Resid（）以及上述单元格的转换。这使我们（再次转移后），带有（混淆）显示的可估计功能。该矩阵有4行，但其等级仅为2行，如该结果的QR分解所确定。然后，我提取R部分的2x2部分，以完成 f 统计量的计算。

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

这里的基本想法是，我需要将所有对比的矩阵分解为两个部分 - 已知的可估计功能和剩余的剩余。一个有趣的方面是，后者的等级是已知，这是我没有利用的事实。在未来的开发中，使用SVD，而不是与qr.isid（）使用SVD，而不是这些旋转。总有新的东西要学习...

原文

I have relied on the qr() function a lot in dealing with rank-deficient situations, but have recently run into some examples where it doesn't work correctly. Consider the
matrix badX below:

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))

We cannot invert this matrix using the solve():

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

Yet qr() and its associated routines thinks this matrix has a rank of 4 and it can invert it:

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

This is a pretty ugly result. I have tried varying the tol argument, with no change in results.

For context, the origin of this result is this contrast matrix:

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")))

... from which I obtained the QR decomposition of its transpose, to find that it is supposedly of rank 4:

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

The above matrix badX is equal to qr.R(badQR)[1:4, 1:4] which based on the rank calculation, was supposed to be a full-rank upper-triangular matrix.

My remedy seems to be to use zapsmall() so that I get the rank right...

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

My question is, why does this happen? If you look at badL, it's pretty clear that it has three zero rows and only the second row is nonzero. I would have thought that qr()'s pivoting methods would work better with this. Is there a better way to get more reliable code?

I am running Windows 11 Pro, version 10.0.22000 build 22000. Here's my R system information.

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"
## 
## I have relied on the qr() function a lot in dealing with rank-deficient situations, but have recently run into some examples where it doesn't work correctly. Consider the

matrix badX below:
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))

We cannot invert this matrix using the solve():
solve(badX)
## Error in solve.default(badX): system is computationally singular: reciprocal condition number = 5.55308e-18

Yet qr() and its associated routines thinks this matrix has a rank of 4 and it can invert it:
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

This is a pretty ugly result. I have tried varying the tol argument, with no change in results.
For context, the origin of this result is this contrast matrix:
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")))

... from which I obtained the QR decomposition of its transpose, to find that it is supposedly of rank 4:
badQR <- qr(t(badL))
badQR$rank
## [1] 4

The above matrix badX is equal to qr.R(badQR)[1:4, 1:4] which based on the rank calculation, was supposed to be a full-rank upper-triangular matrix.
My remedy seems to be to use zapsmall() so that I get the rank right...
qr(zapsmall(t(badL)))$rank
## [1] 1

My question is, why does this happen? If you look at badL, it's pretty clear that it has three zero rows and only the second row is nonzero. I would have thought that qr()'s pivoting methods would work better with this. Is there a better way to get more reliable code?
I am running Windows 11 Pro, version 10.0.22000 build 22000. Here's my R system information.
svn rev`
## [1] "82229"
## 
## $language
## [1] "R"
## 
## $version.string
## [1] "R version 4.2.0 (2022-04-22 ucrt)"
## 
## $nickname
## [1] "Vigorous Calisthenics"

^{Created on 2022-06-21 by the reprex package (v2.0.1)}

More on context

This question came up because I was trying to produce results like this (for a simpler example) in the emmeans package:

> (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

... and in particular, the (confounded) part. This example is with a two-factor model with wool at 2 levels and tension at 3 levels; however, one of the factor combinations is omitted in the data, which means that we can estimate only 1 d.f. for each of the tension main effect and the wool:tension interaction effect, and no main effect at all for wool. There being 4 d.f. for all possible contrasts of the 5 nonempty cells, there are 2 d.f. left over, and those are in the confounded) part.

The computation is based on the associated estimable functions:

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

I have relied on the qr() function a lot in dealing with rank-deficient situations, but have recently run into some examples where it doesn't work correctly. Consider the

matrix badX below:
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))

We cannot invert this matrix using the solve():
solve(badX)
## Error in solve.default(badX): system is computationally singular: reciprocal condition number = 5.55308e-18

Yet qr() and its associated routines thinks this matrix has a rank of 4 and it can invert it:
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

This is a pretty ugly result. I have tried varying the tol argument, with no change in results.
For context, the origin of this result is this contrast matrix:
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")))

... from which I obtained the QR decomposition of its transpose, to find that it is supposedly of rank 4:
badQR <- qr(t(badL))
badQR$rank
## [1] 4

The above matrix badX is equal to qr.R(badQR)[1:4, 1:4] which based on the rank calculation, was supposed to be a full-rank upper-triangular matrix.
My remedy seems to be to use zapsmall() so that I get the rank right...
qr(zapsmall(t(badL)))$rank
## [1] 1

My question is, why does this happen? If you look at badL, it's pretty clear that it has three zero rows and only the second row is nonzero. I would have thought that qr()'s pivoting methods would work better with this. Is there a better way to get more reliable code?
I am running Windows 11 Pro, version 10.0.22000 build 22000. Here's my R system information.
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"
## 
## I have relied on the qr() function a lot in dealing with rank-deficient situations, but have recently run into some examples where it doesn't work correctly. Consider the

matrix badX below:
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))

We cannot invert this matrix using the solve():
solve(badX)
## Error in solve.default(badX): system is computationally singular: reciprocal condition number = 5.55308e-18

Yet qr() and its associated routines thinks this matrix has a rank of 4 and it can invert it:
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

This is a pretty ugly result. I have tried varying the tol argument, with no change in results.
For context, the origin of this result is this contrast matrix:
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")))

... from which I obtained the QR decomposition of its transpose, to find that it is supposedly of rank 4:
badQR <- qr(t(badL))
badQR$rank
## [1] 4

The above matrix badX is equal to qr.R(badQR)[1:4, 1:4] which based on the rank calculation, was supposed to be a full-rank upper-triangular matrix.
My remedy seems to be to use zapsmall() so that I get the rank right...
qr(zapsmall(t(badL)))$rank
## [1] 1

My question is, why does this happen? If you look at badL, it's pretty clear that it has three zero rows and only the second row is nonzero. I would have thought that qr()'s pivoting methods would work better with this. Is there a better way to get more reliable code?
I am running Windows 11 Pro, version 10.0.22000 build 22000. Here's my R system information.
svn rev`
## [1] "82229"
## 
## $language
## [1] "R"
## 
## $version.string
## [1] "R version 4.2.0 (2022-04-22 ucrt)"
## 
## $nickname
## [1] "Vigorous Calisthenics"

^{Created on 2022-06-21 by the reprex package (v2.0.1)}
More on context
This question came up because I was trying to produce results like this (for a simpler example) in the emmeans package:
> (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

... and in particular, the (confounded) part. This example is with a two-factor model with wool at 2 levels and tension at 3 levels; however, one of the factor combinations is omitted in the data, which means that we can estimate only 1 d.f. for each of the tension main effect and the wool:tension interaction effect, and no main effect at all for wool. There being 4 d.f. for all possible contrasts of the 5 nonempty cells, there are 2 d.f. left over, and those are in the confounded) part.
The computation is based on the associated estimable functions:
wool:tension`
     (Intercept) woolB tensionM tensionH woolB:tensionM woolB:tensionH
[1,]           0     0        0        0              1              0

I have relied on the qr() function a lot in dealing with rank-deficient situations, but have recently run into some examples where it doesn't work correctly. Consider the

matrix badX below:
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))

We cannot invert this matrix using the solve():
solve(badX)
## Error in solve.default(badX): system is computationally singular: reciprocal condition number = 5.55308e-18

Yet qr() and its associated routines thinks this matrix has a rank of 4 and it can invert it:
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

This is a pretty ugly result. I have tried varying the tol argument, with no change in results.
For context, the origin of this result is this contrast matrix:
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")))

... from which I obtained the QR decomposition of its transpose, to find that it is supposedly of rank 4:
badQR <- qr(t(badL))
badQR$rank
## [1] 4

The above matrix badX is equal to qr.R(badQR)[1:4, 1:4] which based on the rank calculation, was supposed to be a full-rank upper-triangular matrix.
My remedy seems to be to use zapsmall() so that I get the rank right...
qr(zapsmall(t(badL)))$rank
## [1] 1

My question is, why does this happen? If you look at badL, it's pretty clear that it has three zero rows and only the second row is nonzero. I would have thought that qr()'s pivoting methods would work better with this. Is there a better way to get more reliable code?
I am running Windows 11 Pro, version 10.0.22000 build 22000. Here's my R system information.
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"
## 
## I have relied on the qr() function a lot in dealing with rank-deficient situations, but have recently run into some examples where it doesn't work correctly. Consider the

matrix badX below:
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))

We cannot invert this matrix using the solve():
solve(badX)
## Error in solve.default(badX): system is computationally singular: reciprocal condition number = 5.55308e-18

Yet qr() and its associated routines thinks this matrix has a rank of 4 and it can invert it:
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

This is a pretty ugly result. I have tried varying the tol argument, with no change in results.
For context, the origin of this result is this contrast matrix:
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")))

... from which I obtained the QR decomposition of its transpose, to find that it is supposedly of rank 4:
badQR <- qr(t(badL))
badQR$rank
## [1] 4

The above matrix badX is equal to qr.R(badQR)[1:4, 1:4] which based on the rank calculation, was supposed to be a full-rank upper-triangular matrix.
My remedy seems to be to use zapsmall() so that I get the rank right...
qr(zapsmall(t(badL)))$rank
## [1] 1

My question is, why does this happen? If you look at badL, it's pretty clear that it has three zero rows and only the second row is nonzero. I would have thought that qr()'s pivoting methods would work better with this. Is there a better way to get more reliable code?
I am running Windows 11 Pro, version 10.0.22000 build 22000. Here's my R system information.
svn rev`
## [1] "82229"
## 
## $language
## [1] "R"
## 
## $version.string
## [1] "R version 4.2.0 (2022-04-22 ucrt)"
## 
## $nickname
## [1] "Vigorous Calisthenics"

^{Created on 2022-06-21 by the reprex package (v2.0.1)}
More on context
This question came up because I was trying to produce results like this (for a simpler example) in the emmeans package:
> (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

... and in particular, the (confounded) part. This example is with a two-factor model with wool at 2 levels and tension at 3 levels; however, one of the factor combinations is omitted in the data, which means that we can estimate only 1 d.f. for each of the tension main effect and the wool:tension interaction effect, and no main effect at all for wool. There being 4 d.f. for all possible contrasts of the 5 nonempty cells, there are 2 d.f. left over, and those are in the confounded) part.
The computation is based on the associated estimable functions:
(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

... and on the contrasts among all cells in the design:

> 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

The method I use is to combine the estimable functions for tension and wool:tension, and obtain the QR decomposition of its transpose. Then I use qr.resid() with that and the transpose of the above cell contrasts. That leaves us (after transposing yet again) with the estimable functions shown for (confounded). That matrix has 4 rows, but its rank is only 2, as determined by the QR decomposition of this result; then I extract the 2x2 portion of the R part to complete the computation of the F statistic.

The example at the beginning of this question is similar, but with a larger, more complex model; the badL matrix is the result of the qr.resid() procedure described above. In this context, some of those rows arguably should be zero. My workaround at present is to examine the diagonal elements of R (badR in the OP) and select those that exceed an absolute threshold.

The essential idea here is that I need to decompose that matrix of all contrasts into two parts -- the known estimable functions and the leftovers. And an interesting aspect of this is that the rank of the latter part is known, a fact that I have not taken advantage of. In future development, it may well be, per @duffymo, to use a SVD rather than these gyrations with qr.resid(). There's always new stuff to learn...

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

嗫嚅 2025-02-16 07:06:45

您抱怨说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）根据定义，基于）和svd do。），

x x x

如果我们通过正确乘以完整的级别来重新尺寸 x x x 对角矩阵 d ，QR分解不会改变。

xd = qr 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 Norm，2-Norm）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;

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（大约1E-16）之间的任何阈值是合理的。使用.machine $ double.eps恢复常规QR结果，使您排名4。

“ SQRT”阈值来自我们希望查看X'X的情况。，将x的条件号码平方。

You are complaining that solve can not invert a matrix that seems to be full rank (according to qr). And you think that solve is doing the correct thing, while qr is not.

Well, do not trust solve. It is not a robust numerical procedure and we can fool it easily. Here is a diagonal matrix. It is certainly invertible (by simply inverting its diagonal elements), but solve just can't do it.

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

Now let's look at your badX, which I call R (as it is the upper triangular matrix returned by QR factorization).

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 can not invert it, but qr.solve gives you a proper inverse matrix.

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 factorization is numerically stable, by being less sensitive to the scale (or size, magnitude) of different columns. (Other factorizations like LU (on which solve is based) and SVD do.) By definition, this factorization does

X = Q R

If we re-scale X's columns by right multiplying a full-rank diagonal matrix D, the QR factorization does not change.

X D = Q R D

So let's look at your big matrix t(badL) to which you apply the QR factorization. I call it 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

Let's re-scale its columns so that every column has Euclidean norm (L2 norm, 2-norm) 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

What do you think now? Is it still a matrix with only one nonzero column? Although qr(X) does not actually first re-scale all columns before QR factorization, looking at XD does help you better understand why QR factorization is more robust.

If you do want to intervene, do not use zapsmall; threshold columns by their 2-norm, instead.

X0 <- X
X0[, norm2 < max(norm2) * sqrt(.Machine$double.eps)] <- 0
QR0 <- qr(X0)

QR0$rank
# [1] 1

How do we know that sqrt(.Machine$double.eps) is an appropriate threshold?

Any threshold between sqrt(.Machine$double.eps) (about 1e-8) and .Machine$double.eps (about 1e-16) is reasonable. Using .Machine$double.eps recovers the regular QR result, giving you rank 4.

The "sqrt" threshold comes from the situation where we want to look at X'X, which squares the condition number of X.

回复收藏 0 原文