Compute Approximate CONDition number and 1-Norm of (Large) Matrices

“Estimate”, i.e. compute approximately the CONDition number of a (potentially large, often sparse) matrix A. It works by apply a fast randomized approximation of the 1-norm, norm(A,"1"), through onenormest(.).

Usage

condest(A, t = min(n, 5), normA = norm(A, "1"),
        silent = FALSE, quiet = TRUE)

onenormest(A, t = min(n, 5), A.x, At.x, n,
           silent = FALSE, quiet = silent,
           iter.max = 10, eps = 4 * .Machine$double.eps)

Arguments

A: a square matrix, optional for onenormest(), where instead of A, A.x and At.x can be specified, see there.
t: number of columns to use in the iterations.
normA: number; (an estimate of) the 1-norm of A, by default norm(A, "1"); may be replaced by an estimate.
silent: logical indicating if warning and (by default) convergence messages should be displayed.
quiet: logical indicating if convergence messages should be displayed.
A.x, At.x: when A is missing, these two must be given as functions which compute A %% x, or t(A) %% x, respectively.
n: == nrow(A), only needed when A is not specified.
iter.max: maximal number of iterations for the 1-norm estimator.
eps: the relative change that is deemed irrelevant.

Details

condest() calls lu(A), and subsequently onenormest(A.x = , At.x = ) to compute an approximate norm of the inverse of A, \(A^{-1}\), in a way which keeps using sparse matrices efficiently when A is sparse.

Note that onenormest() uses random vectors and hence both functions' results are random, i.e., depend on the random seed, see, e.g., set.seed().

Value

Both functions return a list; condest() with components,

est: a number \(> 0\), the estimated (1-norm) condition number \(\hat\kappa\); when \(r :=\)rcond(A), \(1/\hat\kappa \approx r\).
v: the maximal \(A x\) column, scaled to norm(v) = 1. Consequently, \(norm(A v) = norm(A) / est\); when est is large, v is an approximate null vector.

The function onenormest() returns a list with components,

est: a number \(> 0\), the estimated norm(A, "1").
v: 0-1 integer vector length n, with an 1 at the index j with maximal column A[,j] in \(A\).
w: numeric vector, the largest \(A x\) found.
iter: the number of iterations used.

References

Nicholas J. Higham and Françoise Tisseur (2000). A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra. SIAM J. Matrix Anal. Appl. 21, 4, 1185–1201.

William W. Hager (1984). Condition Estimates. SIAM J. Sci. Stat. Comput. 5, 311–316.

Author

This is based on octave's condest() and onenormest() implementations with original author Jason Riedy, U Berkeley; translation to R and adaption by Martin Maechler.

Examples

data(KNex, package = "Matrix")
mtm <- with(KNex, crossprod(mm))
system.time(ce <- condest(mtm))
#>    user  system elapsed 
#>   0.043   0.000   0.043 
sum(abs(ce$v)) ## || v ||_1  == 1
#> [1] 1
## Prove that  || A v || = || A || / est  (as ||v|| = 1):
stopifnot(all.equal(norm(mtm %*% ce$v),
                    norm(mtm) / ce$est))

## reciprocal
1 / ce$est
#> [1] 8.123888e-06
system.time(rc <- rcond(mtm)) # takes ca  3 x  longer
#> Warning: 'rcond' via sparse -> dense coercion
#>    user  system elapsed 
#>   0.013   0.006   0.012 
rc
#> [1] 8.123888e-06
all.equal(rc, 1/ce$est) # TRUE -- the approximation was good
#> [1] TRUE

one <- onenormest(mtm)
str(one) ## est = 12.3
#> List of 4
#>  $ est : num 12.3
#>  $ v   : int [1:712] 0 0 0 1 0 0 0 0 0 0 ...
#>  $ w   : num [1:712] 0.036 0 0.0433 0.0375 0 ...
#>  $ iter: int 3
## the maximal column:
which(one$v == 1) # mostly 4, rarely 1, depending on random seed
#> [1] 4