R kappa 计算或估计矩阵的条件数

说明

正则(方)矩阵的条件数是矩阵范数与其逆矩阵范数(或pseudo-inverse)的乘积，因此取决于matrix-norm的类型。

kappa() 默认计算(估计)矩阵或 QR 分解的 R 矩阵的 2 范数条件数(可能是线性拟合)。 2-范数条件数可以表示为矩阵的最大非零奇异值与最小非零奇异值之比。

rcond() 计算条件数倒数的近似值，请参阅详细信息。

用法

kappa(z, ...)
## Default S3 method:
kappa(z, exact = FALSE,
      norm = NULL, method = c("qr", "direct"),
      inv_z = solve(z),
      triangular = FALSE, uplo = "U", ...)

## S3 method for class 'lm'
kappa(z, ...)
## S3 method for class 'qr'
kappa(z, ...)

.kappa_tri(z, exact = FALSE, LINPACK = TRUE, norm = NULL, uplo = "U", ...)

rcond(x, norm = c("O","I","1"), triangular = FALSE, uplo = "U", ...)

参数

z , x	数字或复数矩阵或 qr 的结果或来自从 "lm" 继承的类的拟合。
exact	合乎逻辑的。结果应该是精确的(不超过小舍入误差)而不是快速的(但相当不准确)？
norm	字符串，指定要计算条件数的矩阵范数，请参阅函数 norm() 。对于 kappa() ，默认值为 "2" ，对于 rcond() ，默认值为 "O" ，对于 .kappa_tri() ，默认值取决于 exact ：如果为真，则默认值为 "2" ，否则"O" ，表示一范数或一范数。对于 exact=FALSE ，当前唯一可能的值是无穷范数的 "I" 。对于 exact=TRUE ，范数可以是 "2" ，或 norm(., type = *) 中任何可能的 type 值。
method	部分匹配的字符串，指定要使用的方法； "qr" 主要是向后兼容的默认值。
inv_z	对于 exact=TRUE, norm != "2" ，(近似值)solve(z) ；可以是矩阵 z 的伪逆矩阵或快速近似逆矩阵。默认情况下，当 exact 为 true 时，solve(z) 是条件计算中最昂贵的部分。
triangular	合乎逻辑的。如果为 true，则使用的矩阵只是 z (或 x )的上三角部分或下三角部分，具体取决于
uplo	字符串，"U" 或 "L" 。仅在 triangular = TRUE 时使用，指示是否使用矩阵的上三角部分或下三角部分。
LINPACK	合乎逻辑的。如果为 true 并且 z 不复杂，则调用 LINPACK 例程 dtrco()；否则相关的 LAPACK 例程是。
...	传入或传出其他方法的进一步参数；对于 kappa.*() ，特别是当 norm 不是 "2" 时的 LINPACK 。

细节

对于 kappa() ，如果 exact = FALSE (默认)，则通过 qr(x) 分解 z = QR 的三角矩阵 R 的 1-范数的廉价近似来估计条件数。然而，精确的 2-范数计算(通过 svd )也可能足够快。

请注意，近似 1- 和 Inf-norm 条件数通过method = "direct"计算速度要快得多，并且rcond()计算这些r互惠的条件数，也适用于复杂矩阵，使用标准 LAPACK 例程。目前，还kappa*()函数计算这些近似值exact为 false，即默认情况下。

kappa 和 rcond 是部分相同函数的不同接口。

.kappa_tri是一个内部函数，由kappa.qr和kappa.default;tri是为了三Angular 及其方法仅考虑矩阵的上三角部分或下三角部分，具体取决于uplo = "U"或者"L"，其中"U"之前是内部硬连线的R4.4.0。

底层 LAPACK 代码的不成功结果将导致错误，并给出正错误代码：这些只能通过详细研究 FORTRAN 代码来解释。

值

条件编号 kappa 或近似值 exact = FALSE 。

例子

kappa(x1 <- cbind(1, 1:10)) # 15.71
kappa(x1, exact = TRUE)     # 13.68
kappa(x2 <- cbind(x1, 2:11)) # high! [x2 is singular!]

hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, `+`) }
sv9 <- svd(h9 <- hilbert(9))$ d
kappa(h9)  # pretty high; by default {exact=FALSE, method="qr"} :
kappa(h9) == kappa(qr.R(qr(h9)), norm = "1")
all.equal(kappa(h9, exact = TRUE), # its definition:
          max(sv9) / min(sv9),
          tolerance = 1e-12) ## the same (typically down to 2.22e-16)
kappa(h9, exact = TRUE) / kappa(h9)  # 0.677 (i.e., rel.error = 32%)

## Exact kappa for rectangular matrix
## panmagic.6npm1(7) :
pm7 <- rbind(c( 1, 13, 18, 23, 35, 40, 45),
             c(37, 49,  5, 10, 15, 27, 32),
             c(24, 29, 41, 46,  2, 14, 19),
             c(11, 16, 28, 33, 38, 43,  6),
             c(47,  3,  8, 20, 25, 30, 42),
             c(34, 39, 44,  7, 12, 17, 22),
             c(21, 26, 31, 36, 48,  4,  9))

kappa(pm7, exact=TRUE, norm="1") # no problem for square matrix

m76 <- pm7[,1:6]
(m79 <- cbind(pm7, 50:56, 63:57))

## Moore-Penrose inverse { ~= MASS::ginv(); differing tol (value & meaning)}:
## pinv := p(seudo) inv(erse)
pinv <- function(X, s = svd(X), tol = 64*.Machine$double.eps) {
    if (is.complex(X))
        s$u <- Conj(s$u)
    dx <- dim(X)
    ## X = U D V' ==> Result =  V {1/D} U'
    pI <- function(u,d,v) tcrossprod(v, u / rep(d, each = dx[1L]))
    pos <- (d <- s$d) > max(tol * max(dx) * d[1L], 0)
    if (all(pos))
        pI(s$u, d, s$v)
    else if (!any(pos))
        array(0, dX[2L:1L])
    else { # some pos, some not:
        i <- which(pos)
        pI(s$u[, i, drop = FALSE], d[i],
           s$v[, i, drop = FALSE])
    }
}

## rectangular
kappa(m76, norm="1")
try( kappa(m76, exact=TRUE, norm="1") )# error in  solve().. must be square

## ==> use pseudo-inverse instead of solve() for rectangular {and norm != "2"}:
iZ <- pinv(m76)
kappa(m76, exact=TRUE, norm="1", inv_z = iZ)
kappa(m76, exact=TRUE, norm="M", inv_z = iZ)
kappa(m76, exact=TRUE, norm="I", inv_z = iZ)

iX <- pinv(m79)
kappa(m79, exact=TRUE, norm="1", inv_z = iX)
kappa(m79, exact=TRUE, norm="M", inv_z = iX)
kappa(m79, exact=TRUE, norm="I", inv_z = iX)

作者

The design was inspired by (but differs considerably from) the S function of the same name described in Chambers (1992).

来源

LAPACK 例程DTRCON 和ZTRCON 以及LINPACK 例程DTRCO 。

LAPACK 和 LINPACK 来自https://netlib.org/lapack/ 和https://netlib.org/linpack/，它们的指南在参考文献中列出。

参考

Anderson. E. and ten others (1999) LAPACK Users' Guide. Third Edition. SIAM.
Available on-line at https://netlib.org/lapack/lug/lapack_lug.html.

Chambers, J. M. (1992) Linear models. Chapter 4 of Statistical Models in S eds J. M. Chambers and T. J. Hastie, Wadsworth & Brooks/Cole.

Dongarra, J. J., Bunch, J. R., Moler, C. B. and Stewart, G. W. (1978) LINPACK Users Guide. Philadelphia: SIAM Publications.

也可以看看

norm; svd 用于奇异值分解，qr 用于QR 分解。