R rankMatrix 矩阵的秩

说明

计算‘the’矩阵秩，理论上是一个定义明确的函数(*)，但在实践中有些模糊。我们提供了几种方法，默认对应Matlab的定义。

(*) n \times m 矩阵 A 、 rk(A) 的秩是线性独立列(或行)的最大数量；因此rk(A) \le min(n,m)。

用法

rankMatrix(x, tol = NULL,
           method = c("tolNorm2", "qr.R", "qrLINPACK", "qr",
                      "useGrad", "maybeGrad"),
           sval = svd(x, 0, 0)$d, warn.t = TRUE, warn.qr = TRUE)

qr2rankMatrix(qr, tol = NULL, isBqr = is.qr(qr), do.warn = TRUE)

参数

x	比如说，维度为 n \times m 的数字矩阵。
tol	非负数，指定用于测试 “practically zero” 的(相对，“scalefree”)容差，具体含义取决于 method ；默认情况下， max(dim(x)) * .Machine$double.eps 是根据 Matlab 的默认值(因为它唯一的方法是我们的 method="tolNorm2" )。
method	指定排名计算方法的字符串，可以缩写为： "tolNorm2"： 奇异值的数量>= tol * max(sval)； "qrLINPACK"： 对于稠密矩阵，这是qr(x, tol, LAPACK=FALSE)(这是qr(...)$rank); 这个("qr*"，密集)版本曾经是 过去一段时间推荐的计算矩阵秩的方法。 对于稀疏 x ，这相当于 "qr.R" 。 "qr.R"： 这是三角矩阵 R 的秩，其中 qr() 使用 LAPACK 或 "sparseQR" 方法(参见 qr-methods )来计算分解 QR 。 R 的排名然后定义为 R 的 “non-zero” 对角线条目 d_i 的数量，并且 “non-zero” 满足 |d_i| \ge \mathrm{tol}\cdot\max(|d_i|) 。 "qr"： 用于向后兼容；对于密集 x ，它对应于 "qrLINPACK" ，而对于稀疏 x ，它使用 "qr.R" 。 对于所有 "qr*" 方法，不使用奇异值 sval，这对于大型稀疏矩阵 x 可能至关重要，在这种情况下，当未指定 sval 时，默认计算 svd() 目前将 x 强制转换为密集矩阵。 "useGrad"： 考虑到(递减的)奇异值的“gradient”，最小间隙的索引。 "maybeGrad"： 仅当看起来合理时才选择方法"useGrad"；否则使用"tolNorm2"。
sval	x 的非递增奇异值的数值向量；通常未指定，并在需要时根据 x 计算，即除非 method = "qr" 。
warn.t	逻辑指示 rankMatrix() 在需要 t(x) 而不是 x 时是否应发出警告。目前，仅对于method = "qr"，默认情况下会发出警告，因为调用者通常可以更有效地直接传递t(x)。
warn.qr	在 QR 情况下(即，如果 method 以 "qr" 开头)， rankMatrix() 调用 qr2rankMarix(.., do.warn = warn.qr) ，请参见下文。
qr	一个R对象产生于qr(x,..)，即通常继承自class "qr"或者"sparseQR".
isBqr	logical 指示 qr 是否由 base qr() 产生。 (否则，它通常来自 Matrix 包稀疏 qr 。)
do.warn	逻辑性；如果为 true，则警告 QR 分解的 R 矩阵中的非有限(或在 sparseQR 情况下为负)对角线条目。不要轻易改变！

细节

qr2rankMatrix() 通常从 rankMatrix() 调用 "qr" * method ，但可以直接使用 - 如果 qr 分解可用，则效率更高。

值

如果x是所有0(或零维度)的矩阵，则秩为零；否则，通常是 1:min(dim(x)) 中的正整数，其属性详细说明所使用的方法。

在极少数情况下，稀疏QR分解“fails”在R、d_i(见上文)的对角线条目中以非有限(通常是NaN条目)结束。然后，发出警告信号(除非 warn.qr /do.warn 不为 true)并返回 NA (具体来说， NA_integer_ )。

注意

对于大型稀疏矩阵 x ，除非您可以自己指定 sval ，否则当前 method = "qr" 可能是唯一可行的矩阵，因为其他矩阵需要 sval 并调用 svd() ，当前将 x 强制为 denseMatrix 这可能非常慢或不可能，具体取决于矩阵维度。

请注意，在稀疏 x 、 method = "qr" 的情况下，计算所有非严格零对角线条目 d_i ，直至包括 Matrix 版本 1.1-0，即该方法隐式使用 tol = 0 ，请参阅还有下面的set.seed(42) 示例。

例子

rankMatrix(cbind(1, 0, 1:3)) # 2

(meths <- eval(formals(rankMatrix)$method))

## a "border" case:
H12 <- Hilbert(12)
rankMatrix(H12, tol = 1e-20) # 12;  but  11  with default method & tol.
sapply(meths, function(.m.) rankMatrix(H12, method = .m.))
## tolNorm2   qr.R  qrLINPACK   qr  useGrad maybeGrad
##       11     11         12   12       11        11
## The meaning of 'tol' for method="qrLINPACK" and *dense* x is not entirely "scale free"
rMQL <- function(ex, M) rankMatrix(M, method="qrLINPACK",tol = 10^-ex)
rMQR <- function(ex, M) rankMatrix(M, method="qr.R",     tol = 10^-ex)
sapply(5:15, rMQL, M = H12) # result is platform dependent
##  7  7  8 10 10 11 11 11 12 12 12  {x86_64}
sapply(5:15, rMQL, M = 1000 * H12) # not identical unfortunately
##  7  7  8 10 11 11 12 12 12 12 12
sapply(5:15, rMQR, M = H12)
##  5  6  7  8  8  9  9 10 10 11 11
sapply(5:15, rMQR, M = 1000 * H12) # the *same*


## "sparse" case:
M15 <- kronecker(diag(x=c(100,1,10)), Hilbert(5))
sapply(meths, function(.m.) rankMatrix(M15, method = .m.))
#--> all 15, but 'useGrad' has 14.
sapply(meths, function(.m.) rankMatrix(M15, method = .m., tol = 1e-7)) # all 14

## "large" sparse
n <- 250000; p <- 33; nnz <- 10000
L <- sparseMatrix(i = sample.int(n, nnz, replace=TRUE),
                  j = sample.int(p, nnz, replace=TRUE),
                  x = rnorm(nnz))
(st1 <- system.time(r1 <- rankMatrix(L)))                # warning+ ~1.5 sec (2013)
(st2 <- system.time(r2 <- rankMatrix(L, method = "qr"))) # considerably faster!
r1[[1]] == print(r2[[1]]) ## -->  ( 33  TRUE )

## another sparse-"qr" one, which ``failed'' till 2013-11-23:
set.seed(42)
f1 <- factor(sample(50, 1000, replace=TRUE))
f2 <- factor(sample(50, 1000, replace=TRUE))
f3 <- factor(sample(50, 1000, replace=TRUE))
D <- t(do.call(rbind, lapply(list(f1,f2,f3), as, 'sparseMatrix')))
dim(D); nnzero(D) ## 1000 x 150 // 3000 non-zeros (= 2%)
stopifnot(rankMatrix(D,           method='qr') == 148,
	  rankMatrix(crossprod(D),method='qr') == 148)

## zero matrix has rank 0 :
stopifnot(sapply(meths, function(.m.)
                        rankMatrix(matrix(0, 2, 2), method = .m.)) == 0)

作者

Martin Maechler; for the "*Grad" methods building on suggestions by Ravi Varadhan.

也可以看看

qr、svd。