Estimate the Reciprocal Condition Number

Estimate the reciprocal of the condition number of a matrix.

This is a generic function with several methods, as seen by showMethods(rcond).

Usage

rcond(x, norm, ...)

# S4 method for class 'sparseMatrix,character'
rcond(x, norm, useInv=FALSE, ...)

Arguments

x

an R object that inherits from the Matrix class.

norm

character string indicating the type of norm to be used in the estimate. The default is "O" for the 1-norm ("O" is equivalent to "1"). For sparse matrices, when useInv=TRUE, norm can be any of the kinds allowed for norm; otherwise, the other possible value is "I" for the infinity norm, see also norm.

useInv

logical (or "Matrix" containing solve(x)). If not false, compute the reciprocal condition number as $1/(\|x\| \cdot \|x^{-1}\|)$, where $x^{-1}$ is the inverse of $x$, solve(x).

This may be an efficient alternative (only) in situations where solve(x) is fast (or known), e.g., for (very) sparse or triangular matrices.

Note that the result may differ depending on useInv, as per default, when it is false, an approximation is computed.

...

further arguments passed to or from other methods.

Value

An estimate of the reciprocal condition number of x.

BACKGROUND

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

More generally, the condition number is defined (also for non-square matrices $A$) as $$\kappa(A) = \frac{\max_{\|v\| = 1} \|A v\|}{\min_{\|v\| = 1} \|A v\|}.$$ Whenever x is not a square matrix, in our method definitions, this is typically computed via rcond(qr.R(qr(X)), ...) where X is x or t(x).

The condition number takes on values between 1 and infinity, inclusive, and can be viewed as a factor by which errors in solving linear systems with this matrix as coefficient matrix could be magnified.

rcond() computes the reciprocal condition number $1/\kappa$ with values in $[0,1]$ and can be viewed as a scaled measure of how close a matrix is to being rank deficient (aka “singular”).

Condition numbers are usually estimated, since exact computation is costly in terms of floating-point operations. An (over) estimate of reciprocal condition number is given, since by doing so overflow is avoided. Matrices are well-conditioned if the reciprocal condition number is near 1 and ill-conditioned if it is near zero.

References

Golub, G., and Van Loan, C. F. (1989). Matrix Computations, 2nd edition, Johns Hopkins, Baltimore.

Examples

x <- Matrix(rnorm(9), 3, 3)
rcond(x)
#> [1] 0.02377467
## typically "the same" (with more computational effort):
1 / (norm(x) * norm(solve(x)))
#> [1] 0.02377467
rcond(Hilbert(9))  # should be about 9.1e-13
#> [1] 9.093806e-13

## For non-square matrices:
rcond(x1 <- cbind(1,1:10))# 0.05278
#> [1] 0.05278005
rcond(x2 <- cbind(x1, 2:11))# practically 0, since x2 does not have full rank
#> [1] 1.144414e-17

## sparse
(S1 <- Matrix(rbind(0:1,0, diag(3:-2))))
#> 8 x 6 sparse Matrix of class "dgCMatrix"
#>                   
#> [1,] . 1 . 1  .  1
#> [2,] . . . .  .  .
#> [3,] 3 . . .  .  .
#> [4,] . 2 . .  .  .
#> [5,] . . 1 .  .  .
#> [6,] . . . .  .  .
#> [7,] . . . . -1  .
#> [8,] . . . .  . -2
rcond(S1)
#> Warning: 'rcond' via sparse -> dense coercion
#> [1] 0.2230509
m1 <- as(S1, "denseMatrix")
all.equal(rcond(S1), rcond(m1))
#> Warning: 'rcond' via sparse -> dense coercion
#> [1] TRUE

## wide and sparse
rcond(Matrix(cbind(0, diag(2:-1))))
#> Warning: matrix is structurally rank deficient; using augmented matrix with additional 1 row(s) of zeros
#> Warning: 'rcond' via sparse -> dense coercion
#> [1] 0

## Large sparse example ----------
m <- Matrix(c(3,0:2), 2,2)
M <- bdiag(kronecker(Diagonal(2), m), kronecker(m,m))
36*(iM <- solve(M)) # still sparse
#> 8 x 8 sparse Matrix of class "dtCMatrix"
#>                            
#> [1,] 12 -6  .  . .  .  .  .
#> [2,]  . 18  .  . .  .  .  .
#> [3,]  .  . 12 -6 .  .  .  .
#> [4,]  .  .  . 18 .  .  .  .
#> [5,]  .  .  .  . 4 -2 -2  1
#> [6,]  .  .  .  . .  6  . -3
#> [7,]  .  .  .  . .  .  6 -3
#> [8,]  .  .  .  . .  .  .  9
MM <- kronecker(Diagonal(10), kronecker(Diagonal(5),kronecker(m,M)))
dim(M3 <- kronecker(bdiag(M,M),MM)) # 12'800 ^ 2
#> [1] 12800 12800
if(interactive()) ## takes about 2 seconds if you have >= 8 GB RAM
  system.time(r <- rcond(M3))
## whereas this is *fast* even though it computes  solve(M3)
system.time(r. <- rcond(M3, useInv=TRUE))
#>    user  system elapsed 
#>   0.006   0.000   0.006 
if(interactive()) ## the values are not the same
  c(r, r.)  # 0.05555 0.013888
## for all 4 norms available for sparseMatrix :
cbind(rr <- sapply(c("1","I","F","M"),
             function(N) rcond(M3, norm=N, useInv=TRUE)))
#>           [,1]
#> 1 1.388889e-02
#> I 7.812500e-03
#> F 2.059462e-05
#> M 3.292181e-02