如何根据R中的第二个稀疏矩阵修改大型稀疏矩阵的非零元素

我的方法中是否存在根本错误？

对给你.

# This takes hours.
m2[position_matrix_from_m1[1,], position_matrix_from_m1[2,]] <- 1
m1[position_matrix_from_m1[1,], position_matrix_from_m1[2,]] <- 0

语法as mat[rn, cn](mat是密集矩阵还是稀疏矩阵) Select rn中的所有行和cn中的所有列.得到一个length(rn) x length(cn)的矩阵.下面是一个小例子:

A <- matrix(1:9, 3, 3)
#     [,1] [,2] [,3]
#[1,]    1    4    7
#[2,]    2    5    8
#[3,]    3    6    9
rn <- 1:2
cn <- 2:3
A[rn, cn]
#     [,1] [,2]
#[1,]    4    7
#[2,]    5    8

你要做的是 Select (rc[1], cn[1])(rc[2], cn[2])...，只有正确的语法是mat[cbind(rn, cn)].下面是一个演示:

A[cbind(rn, cn)]
#[1] 4 8

因此，您需要将代码修复为:

m2[cbind(position_matrix_from_m1[1,], position_matrix_from_m1[2,])] <- 1
m1[cbind(position_matrix_from_m1[1,], position_matrix_from_m1[2,])] <- 0

哦，等等...根据你position_matrix_from_m1的构造，这只是

ij <- t(position_matrix_from_m1[1:2, ])
m2[ij] <- 1
m1[ij] <- 0

Now, let me explain how you can do better.你没有充分利用summary().它返回一个3列数据帧，给出(i，j，x)三元组，其中i和j都是从1开始的索引.您可以直接使用这个漂亮的输出，如下所示:

# Getting (i, j, x) triplet (stored as a data.frame) for both `m1` and `m2`
position_matrix_from_m1 <- summary(m1)
# you never seem to use `position_matrix_from_m2` so I skip it

# Subsetting to get the elements of interest.
position_matrix_from_m1 <- subset(position_matrix_from_m1, x > 0 & x < 0.05)

现在你可以做:

ij <- as.matrix(position_matrix_from_m1[, 1:2])
m2[ij] <- 1
m1[ij] <- 0

Is there a even better solution? Yes!注意，m1和m2中的非零元素位于相同的位置.所以基本上，你只需要根据m1@x改变m2@x.

ind <- m1@x > 0 & m1@x < 0.05
m2@x[ind] <- 1
m1@x[ind] <- 0

A complete R session

我没有足够的内存来创建您的大型矩阵，所以我将您的问题大小减少了一点，以便进行测试.一切顺利.

library(Matrix)
# Generating the example matrices.
set.seed(42)

## reduce problem size to what my laptop can bear with
squeeze <- 0.1

# Rows with values.
i <- sample(1:(41000 * squeeze), 227000000 * squeeze ^ 2, replace = TRUE)
# Columns with values.
j <- sample(1:(55000 * squeeze), 227000000 * squeeze ^ 2, replace = TRUE)

# Values for the first matrix.
x1 <- runif(227000000 * squeeze ^ 2)
# Values for the second matrix.
x2 <- sample(1:3, 227000000 * squeeze ^ 2, replace = TRUE)

# Constructing the matrices.
m1 <- sparseMatrix(i = i, j = j, x = x1)
m2 <- sparseMatrix(i = i, j = j, x = x2)

## give me more usable RAM
rm(i, j, x1, x2)

##
## fix to your code
##

m1a <- m1
m2a <- m2

# Getting (i, j, x) triplet (stored as a data.frame) for both `m1` and `m2`
position_matrix_from_m1 <- summary(m1)

# Subsetting to get the elements of interest.
position_matrix_from_m1 <- subset(position_matrix_from_m1, x > 0 & x < 0.05)

ij <- as.matrix(position_matrix_from_m1[, 1:2])
m2a[ij] <- 1
m1a[ij] <- 0

##
## the best solution
##

m1b <- m1
m2b <- m2

ind <- m1@x > 0 & m1@x < 0.05
m2b@x[ind] <- 1
m1b@x[ind] <- 0

##
## they are identical
##
all.equal(m1a, m1b)
#[1] TRUE
all.equal(m2a, m2b)
#[1] TRUE

Caveat:

我知道有些人可能会提议

m1c <- m1
m2c <- m2

logi <- m1 > 0 & m1 < 0.05
m2c[logi] <- 1
m1c[logi] <- 0

它在R的语法中看起来完全自然.但是相信我，对于大型矩阵来说，这是非常缓慢的.