kappa {base} | R Documentation |
The condition number of a regular (square) matrix is the product of the norm of the matrix and the norm of its inverse (or pseudo-inverse), and hence depends on the kind of matrix-norm.
kappa()
computes by default (an estimate of) the 2-norm
condition number of a matrix or of the R
matrix of a QR
decomposition, perhaps of a linear fit. The 2-norm condition number
can be shown to be the ratio of the largest to the smallest
non-zero singular value of the matrix.
rcond()
computes an approximation of the reciprocal
condition number, see the details.
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 |
a numeric or complex matrix or a result of
|
exact |
logical. Should the result be exact (up to small rounding error) as opposed to fast (but quite inaccurate)? |
norm |
character string, specifying the matrix norm with respect
to which the condition number is to be computed, see the function
|
method |
a partially matched character string specifying the method to be used;
|
inv_z |
for |
triangular |
logical. If true, the matrix used is just the upper or
lower triangular part of |
uplo |
character string, either |
LINPACK |
logical. If true and |
... |
further arguments passed to or from other methods;
for |
For kappa()
, if exact = FALSE
(the default) the
condition number is estimated by a cheap approximation to the 1-norm of
the triangular matrix R
of the qr(x)
decomposition
z = QR
. However, the exact 2-norm calculation (via
svd
) is also likely to be quick enough.
Note that the approximate 1- and Inf-norm condition numbers via
method = "direct"
are much faster to
calculate, and rcond()
computes these reciprocal
condition numbers, also for complex matrices, using standard LAPACK
routines.
Currently, also the kappa*()
functions compute these
approximations whenever exact
is false, i.e., by default.
kappa
and rcond
are different interfaces to
partly identical functionality.
.kappa_tri
is an internal function called by kappa.qr
and
kappa.default
; tri
is for triangular and its methods
only consider the upper or lower triangular part of the matrix, depending
on uplo = "U"
or "L"
, where "U"
was internally hard
wired before R 4.4.0.
Unsuccessful results from the underlying LAPACK code will result in an error giving a positive error code: these can only be interpreted by detailed study of the FORTRAN code.
The condition number, kappa
, or an approximation if
exact = FALSE
.
The design was inspired by (but differs considerably from) the S function of the same name described in Chambers (1992).
The LAPACK routines DTRCON
and ZTRCON
and the LINPACK
routine DTRCO
.
LAPACK and LINPACK are from https://netlib.org/lapack/ and https://netlib.org/linpack/ and their guides are listed in the references.
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
for the singular value decomposition and
qr
for the QR
one.
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)