# SCM Repository

[matrix] View of /pkg/Matrix/R/rankMatrix.R
 [matrix] / pkg / Matrix / R / rankMatrix.R

# View of /pkg/Matrix/R/rankMatrix.R

Sun Nov 24 16:02:10 2013 UTC (5 years, 8 months ago) by mmaechler
File size: 4343 byte(s)
`rankMatrix(); method="qr" as it is not related to LINPACK for the important case of sparse x:  It now uses *tol* to count non-zero entries in diag(R)`
```#### Determine *the* rank of a matrix
#### --------------------------------
##
## As this is not such a well-defined problem as people think,
## we provide *some* possibilities here, including the Matlab one.
##
## Ideas by Martin Maechler (April 2007) and Ravi Varadhan (October 2007)

## (for now: use code & examples from /u/maechler/R/MM/NUMERICS/rankMat.R )

rankMatrix <- function(x, tol = NULL,
sval = svd(x, 0,0)\$d, warn.t = TRUE)
{
## Purpose: rank of a matrix ``as Matlab'' or "according to Ravi V"
## ----------------------------------------------------------------------
## Arguments: x: a numerical matrix, maybe non-square
##          tol: numerical tolerance (compared to singular values)
##         sval: vector of non-increasing singular values of  x
##               (pass as argument if already known)
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 7 Apr 2007, 16:16
## ----------------------------------------------------------------------
##
## maybeGrad (Ravi V.): This algorithm determines the rank based on the
## absolute, singular values, rather than enforcing a rigid
## tolerance criterion,
##
## Author: Ravi Varadhan, Date: 22 October 2007 // Tweaks: MM, Oct.23

## ----------------------------------------------------------------------

stopifnot(length(d <- dim(x)) == 2)
p <- min(d)
miss.meth <- missing(method)
## "qrLINPACK" is allowed for backcompatibility [change @ 2013-11-24]
method <- if(!miss.meth && identical(method, "qrLINPACK")) "qr" else match.arg(method)

stopifnot(length(sval) == p,
diff(sval) <= 0) # must be sorted non-increasingly: max = s..[1]
ln.av <- log(abs(sval))
diff1 <- diff(ln.av)
grad <- (min(ln.av) - max(ln.av)) / p
}#  -------
}
x.dense <- is.numeric(x) || is(x,"denseMatrix")
findTol <- function() {         # the "Matlab" default:
stopifnot(diff(sval) <= 0)
max(d) * .Machine\$double.eps * abs(sval[1])
}
if(method == "qr") {
if(is.null(tol)) tol <- findTol() # also for sparse QR
} else { ## (method != "qr")
if(is.null(tol)) {
if(!x.dense && missing(sval) && prod(d) >= 100000L)
warning(gettextf("rankMatrix(<large sparse Matrix>, method = '%s') coerces to dense matrix.\n  Probably should rather use  method = 'qr' !?",
method),
immediate.=TRUE, domain=NA)
tol <- findTol()
} else stopifnot((tol <- as.numeric(tol)[[1]]) >= 0)
}
}

structure(## rank :
else if(method == "qr") {
if(do.t <- (d[1L] < d[2L]))
warning(gettextf(
"rankMatrix(x, method='qr'): computing t(x) as nrow(x) < ncol(x)"))
q.r <- qr(if(do.t) t(x) else x, tol=tol, LAPACK = FALSE)
if(is(q.r, "qr"))
q.r\$rank
else if(is(q.r,"sparseQR")) { # 'tol' not used in qr()
## FIXME: Here, we could be quite a bit faster,
## by not returning the full sparseQR, but just
## doing the following in C, and return the rank.
d <- sort(abs(diag(q.r@R))) ## is abs(.) unneeded? [FIXME]
## declare those entries to be zero that are < tol*max(.)
sum(d > tol * d[length(d)])
## was sum(diag(q.r@R) != 0)
}
else stop(gettextf(
"method %s not applicable for qr() result class %s",
sQuote(method), dQuote(class(q.r)[1])),
domain=NA)
} else sum(sval >= tol),
"method" = method,
"tol" = if(useGrad) NA else tol)
}

## Ravi's plot of the absolute singular values:
if(FALSE) {
## if (plot.eigen) {
plot(abs(sval), type = "b", xlab = "Index", xaxt = "n",
log = "y", ylab = "|singular value|   [log scaled]")
axis(1, at = unique(c(axTicks(1), rank, p)))
abline(v = rank, lty = 3)
mtext(sprintf("rank = %d (used %s (%g))", rank,