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

Description

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

Value

Both functions return a list;

condest() with components,

est: a number $> 0$ , the estimated (1-norm) condition number $\hat{κ}$ ; when $r :=$ rcond(A), $1 / \hat{κ} \approx r$ .
v: the maximal $A x$ column, scaled to norm(v) = 1. Consequently, $n o r m (A v) = n o r m (A) / e s t$ ; 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.

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.

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.

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

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.

Examples

Run this code

data(KNex)
mtm <- with(KNex, crossprod(mm))
system.time(ce <- condest(mtm))
sum(abs(ce$v)) ## || v ||_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
system.time(rc <- rcond(mtm)) # takes ca  3 x  longer
rc
all.equal(rc, 1/ce$est) # TRUE -- the approxmation was good

one <- onenormest(mtm)
str(one) ## est = 12.3
## the maximal column:
which(one$v == 1) # mostly 4, rarely 1, depending on random seed

Run the code above in your browser using DataLab

Last chance! 50% off unlimited learning