1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
|
#### 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)
qr2rankMatrix <- function(qr, tol = NULL, isBqr = is.qr(qr), do.warn=TRUE) {
## NB: 1) base::qr(*, LAPACK = TRUE/FALSE) differ via attr(.,"useLAPACK")
## 2) if LAPACK=TRUE, .$rank is useless (always = full rank)
##
## return ( . ) :
if(isBqr && !isTRUE(attr(qr, "useLAPACK")))
qr$rank
else {
diagR <- if(isBqr) # hence "useLAPACK" here
diag(qr$qr) # faster than, but equivalent to diag(qr.R(q.r))
else ## ==> assume Matrix::qr() i.e., currently "sparseQR"
## 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.
diag(qr@R)
if(anyNA(diagR) || !all(is.finite(diagR))) {
if(do.warn) {
ifi <- is.finite(diagR)
warning(gettextf(
"qr2rankMatrix(.): QR with only %d out of %d finite diag(R) entries",
sum(ifi), length(ifi)))
}
## return
NA_integer_
## alternative: gives *too* small rank in typical cases
## reduce the maximal rank by omitting all non-finite entries:
## diagR <- diagR[is.finite(diagR)]
## if(length(diagR) == 0)
## return(NA_integer_)
} else {
diagR <- abs(diagR) # sign( diag(R) ) are *not* coerced to positive
## declare those entries to be zero that are < tol*max(.)
if((mdi <- max(diagR, na.rm=TRUE)) > 0) {
if(!is.numeric(tol)) {
## d := dim(x) extracted from qr, in both (dense and sparse) qr() cases
d <- dim(if(isBqr) qr$qr else qr)
tol <- max(d) * .Machine$double.eps
}
sum(diagR >= tol * mdi)
## was sum(diag(q.r@R) != 0)
}
else 0L # for 0-matrix or all NaN or negative diagR[]
}
} ## else {Lapack or sparseQR}
}
rankMatrix <- function(x, tol = NULL,
method = c("tolNorm2", "qr.R", "qrLINPACK", "qr",
"useGrad", "maybeGrad"),
sval = svd(x, 0,0)$d, warn.t = TRUE, warn.qr = 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
## "gradient" of 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)
method <- match.arg(method)
if(useGrad <- (method %in% c("useGrad", "maybeGrad"))) {
stopifnot(length(sval) == p)
if(p > 1) stopifnot(diff(sval) <= 0) # must be sorted non-increasingly: max = s..[1]
if(sval[1] == 0) { ## <==> all singular values are zero <==> Matrix = 0 <==> rank = 0
useGrad <- FALSE
method <- eval(formals()[["method"]])[[1]]
} else {
ln.av <- log(abs(sval))
if(any(s0 <- sval == 0)) ln.av[s0] <- - .Machine$double.xmax # so we get diff() == 0
diff1 <- diff(ln.av)
if(method == "maybeGrad") {
grad <- (min(ln.av) - max(ln.av)) / p
useGrad <- !is.na(grad) && p > 1 && min(diff1) <= min(-3, 10 * grad)
}# -------
}
}
if(!useGrad) {
x.dense <- is.numeric(x) || is(x,"denseMatrix")
## "qr" is allowed for backcompatibility [change @ 2013-11-24]
if((Meth <- method) == "qr")
method <- if(x.dense) "qrLINPACK" else "qr.R"
else Meth <- substr(method, 1,2)
if(Meth == "qr") {
if(is.null(tol)) tol <- max(d) * .Machine$double.eps
} else { ## (Meth != "qr"), i.e. "tolNorm2"
if(is.null(tol)) {
if(!x.dense && missing(sval) && prod(d) >= 100000L)
warning(gettextf(
"rankMatrix(<large sparse Matrix>, method = '%s') coerces to dense matrix.
Probably should rather use method = 'qr' !?",
method),
immediate.=TRUE, domain=NA)
## the "Matlab" default:
if(p > 1) stopifnot(diff(sval) <= 0) #=> sval[1]= max(sval)
tol <- max(d) * .Machine$double.eps
} else stopifnot((tol <- as.numeric(tol)[[1]]) >= 0)
}
}
structure(## rank :
if(useGrad) which.min(diff1)
else if(Meth == "qr") {
if((do.t <- (d[1L] < d[2L])) && warn.t)
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 = method != "qrLINPACK")
qr2rankMatrix(q.r, tol=tol, isBqr = x.dense, do.warn = warn.qr)
}
else if(sval[1] > 0) sum(sval >= tol * sval[1]) else 0L, ## "tolNorm2"
"method" = method,
"useGrad" = useGrad,
"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,
if(use.grad)"'gradient'" else "fixed tol.",
if(use.grad) min(diff1) else tol))
}
|
