R程序中非常大的矩阵的svd

Question

我在一个txt文件中有一个矩阵60 000 x 60 000,我需要得到这个矩阵的svd.我使用R但我不知道R是否可以生成它.

Answer 1

我认为svd使用irlba包计算(部分)bigmemory并且bigalgebra不使用大量内存是可能的.

首先让我们创建一个20000*20000矩阵并将其保存到文件中

require(bigmemory)
require(bigalgebra)
require(irlba)

con <- file("mat.txt", open = "a")
replicate(20, {
    x <- matrix(rnorm(1000 * 20000), nrow = 1000)
    write.table(x, file  = 'mat.txt', append = TRUE,
            row.names = FALSE, col.names = FALSE)
})

file.info("mat.txt")$size
## [1] 7.264e+09   7.3 Gb
close(con)

然后你可以使用读取这个矩阵 bigmemory::read.big.matrix

bigm <- read.big.matrix("mat.txt", sep = " ",
                        type = "double",
                        backingfile = "mat.bk",
                        backingpath = "/tmp",
                        descriptorfile = "mat.desc")

str(bigm)
## Formal class 'big.matrix' [package "bigmemory"] with 1 slots
##   ..@ address:<externalptr>

dim(bigm)
## [1] 20000 20000

bigm[1:3, 1:3]
##            [,1]     [,2]     [,3]
## [1,] -0.3623255 -0.58463 -0.23172
## [2,] -0.0011427  0.62771  0.73589
## [3,] -0.1440494 -0.59673 -1.66319

现在我们可以使用优秀的irlba包,如包装插图中所述.

第一步包括定义矩阵乘法运算符,它可以与big.matrix对象一起使用,然后使用该irlba::irlba函数

### vignette("irlba", package = "irlba") # for more info

matmul <- function(A, B, transpose=FALSE) {
    ## Bigalgebra requires matrix/vector arguments
    if(is.null(dim(B))) B <- cbind(B)

    if(transpose)
        return(cbind((t(B) %*% A)[]))

    cbind((A %*% B)[])
}

dim(bigm)

system.time(
S <- irlba(bigm, nu = 2, nv = 2, matmul = matmul)
)

##    user  system elapsed 
## 169.820   0.923 170.194


str(S)
## List of 5
##  $ d    : num [1:2] 283 283
##  $ u    : num [1:20000, 1:2] -0.00615 -0.00753 -0.00301 -0.00615 0.00734 ...
##  $ v    : num [1:20000, 1:2] 0.020086 0.012503 0.001065 -0.000607 -0.006009 ...
##  $ iter : num 10
##  $ mprod: num 310

我忘了设置种子使它可以再现,但我只想表明它可以在R中做到这一点.

编辑

如果您使用的是该软件包的新版本irlba,则上述代码会抛出错误,因为matmult该函数的参数irlba已重命名为mult.因此,您应该更改此部分代码

S <- irlba(bigm, nu = 2, nv = 2, matmul = matmul)

通过

S <- irlba(bigm, nu = 2, nv = 2, mult = matmul)

我要感谢@FrankD指出这一点.

Answer 2

在 R 3.x+ 中，您可以构造一个该大小的矩阵，向量大小的上限为 2^53 （或者可能是 2^53-1 ），高于之前的 2^31-1 ，这就是 Andrie 的原因在他过时的安装上抛出错误。每个数字元素通常需要 10 个字节。好歹：

> 2^53 < 10*60000^2
[1] FALSE  # so you are safe on that account.

它还适合 64GB（但不适合 32GB）：

> 64000000000 < 10*60000^2
[1] FALSE

一般来说，要完成任何严肃的工作，您至少需要最大对象大小的 3 倍，因此即使使用新的扩展向量/矩阵，这似乎也相当临界。