“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 1norm,
norm(A,"1")
, through onenormest(.)
.
1 2 3 4 5 6 
A 
a square matrix, optional for 
t 
number of columns to use in the iterations. 
normA 
number; (an estimate of) the 1norm of 
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 
n 

iter.max 
maximal number of iterations for the 1norm estimator. 
eps 
the relative change that is deemed irrelevant. 
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()
.
Both functions return a list
;
condest()
with components,
est 
a number > 0, the estimated (1norm) condition number
k.; when r := 
v 
the maximal A x column, scaled to norm(v) = 1.
Consequently, norm(A v) = norm(A) / est;
when 
The function onenormest()
returns a list with components,
est 
a number > 0, the estimated 
v 
01 integer vector length 
w 
numeric vector, the largest A x found. 
iter 
the number of iterations used. 
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.
Nicholas J. Higham and Françoise Tisseur (2000). A Block Algorithm for Matrix 1Norm Estimation, with an Application to 1Norm Pseudospectra. SIAM J. Matrix Anal. Appl. 21, 4, 1185–1201. http://dx.doi.org/10.1137/S0895479899356080
William W. Hager (1984). Condition Estimates. SIAM J. Sci. Stat. Comput. 5, 311–316.
norm
, rcond
.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18  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

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
All documentation is copyright its authors; we didn't write any of that.