如何在 R 中的稀疏矩阵中查找和命名连续的非零项?
How to find and name contiguous non-zero entries in a sparse matrix in R?
我的问题在概念上很简单。
我正在寻找一种计算效率高的解决方案(我自己的解决方案附在最后)。
假设我们有一个潜在的非常大的稀疏矩阵,就像下面左边的那个,并且想要 'name' 每个区域的连续非零元素都有一个单独的代码(见右边的矩阵)
1 1 1 . . . . . 1 1 1 . . . . .
1 1 1 . 1 1 . . 1 1 1 . 4 4 . .
1 1 1 . 1 1 . . 1 1 1 . 4 4 . .
. . . . 1 1 . . ---> . . . . 4 4 . .
. . 1 1 . . 1 1 . . 3 3 . . 7 7
1 . 1 1 . . 1 1 2 . 3 3 . . 7 7
1 . . . 1 . . . 2 . . . 5 . . .
1 . . . . 1 1 1 2 . . . . 6 6 6
在我的应用程序中,连续元素将形成矩形、直线或单个点,它们只能通过顶点相互接触(即矩阵中没有 irregular/non 矩形区域)。
我想象的解决方案是将稀疏矩阵表示的行和列索引与具有适当值的向量('name'代码)匹配。我的解决方案使用多个 for loops
并且适用于中小型矩阵,但随着矩阵的维度变大(> 1000),很快就会陷入循环。这可能取决于我在 R 编程方面不是那么先进的事实 - 我找不到任何计算 trick/function 来更好地解决它。
任何人都可以建议一种在 R 中计算更有效的方法吗?
我的解决方案:
mySolution <- function(X){
if (class(X) != "ngCMatrix") {stop("Input must be a Sparse Matrix")}
ind <- which(X == TRUE, arr.ind = TRUE)
r <- ind[,1]
c <- ind[,2]
lr <- nrow(ind)
for (i in 1:lr) {
if(i == 1) {bk <- 1}
else {
if (r[i]-r[i-1] == 1){bk <- c(bk, bk[i-1])}
else {bk <- c(bk, bk[i-1]+1)}
}
}
for (LOOP in 1:(lr-1)) {
tr <- r[LOOP]
tc <- c[LOOP]
for (j in (LOOP+1):lr){
if (r[j] == tr) {
if(c[j] == tc + 1) {bk[j] <- bk[LOOP]}
}
}
}
val <- unique(bk)
for (k in 1:lr){
bk[k] <- which(val==bk[k])
}
return(sparseMatrix(i = r, j = c, x = bk))
}
在此先感谢您的帮助或指点。
在很大程度上依赖于所有要分组的相邻元素仅形成 rectangles/lines/points 这一事实,我们看到矩阵的元素可以根据它们在矩阵上的 [row, col]
索引通过关系聚合(abs(row1 - row2) + abs(col1 - col2)) < 2
。
因此,从 [row, col]
索引开始:
sm = as.matrix(summary(m))
我们计算它们的距离,正如 GiuGe 所说-实际上是 "manhattan" 方法:
d = dist(sm, "manhattan")
Single-linkage 的 属性 在其最近邻元素上的聚类元素在这里很有用。此外,我们可以通过 cutree
ing on "h = 1"(其中索引的距离为“<2”)获得一组元素:
gr = cutree(hclust(d, "single"), h = 1)
最后,我们可以将上面的内容包装在一个新的稀疏矩阵中:
sparseMatrix(i = sm[, "i"], j = sm[, "j"], x = gr)
#8 x 8 sparse Matrix of class "dgCMatrix"
#
#[1,] 1 1 1 . . . . .
#[2,] 1 1 1 . 4 4 . .
#[3,] 1 1 1 . 4 4 . .
#[4,] . . . . 4 4 . .
#[5,] . . 3 3 . . 7 7
#[6,] 2 . 3 3 . . 7 7
#[7,] 2 . . . 5 . . .
#[8,] 2 . . . . 6 6 6
"m"使用的是:
library(Matrix)
m = new("ngCMatrix"
, i = c(0L, 1L, 2L, 5L, 6L, 7L, 0L, 1L, 2L, 0L, 1L, 2L, 4L, 5L, 4L,
5L, 1L, 2L, 3L, 6L, 1L, 2L, 3L, 7L, 4L, 5L, 7L, 4L, 5L, 7L)
, p = c(0L, 6L, 9L, 14L, 16L, 20L, 24L, 27L, 30L)
, Dim = c(8L, 8L)
, Dimnames = list(NULL, NULL)
, factors = list()
)
编辑 2017 年 2 月 10 日
另一个想法(并且再次考虑到相邻元素仅形成 rectangles/lines/points 的事实)是迭代 - 在升序列中 - 通过 [row, col]
索引,并在每一步中找到当前列和行中其最近邻居的每个元素的距离。如果找到 "< 2" 距离,则该元素与其相邻元素分组,否则开始一个新组。包装成一个函数:
ff = function(x)
{
sm = as.matrix(summary(x))
gr = integer(nrow(sm)); ngr = 0L ; gr[1] = ngr
lastSeenRow = integer(nrow(x))
lastSeenCol = integer(ncol(x))
for(k in 1:nrow(sm)) {
kr = sm[k, 1]; kc = sm[k, 2]
i = lastSeenRow[kr]
j = lastSeenCol[kc]
if(i && (abs(kc - sm[i, 2]) == 1)) gr[k] = gr[i]
else if(j && (abs(kr - sm[j, 1]) == 1)) gr[k] = gr[j]
else { ngr = ngr + 1L; gr[k] = ngr }
lastSeenRow[kr] = k
lastSeenCol[kc] = k
}
sparseMatrix(i = sm[, "i"], j = sm[, "j"], x = gr)
}
并应用于 "m":
ff(m)
#8 x 8 sparse Matrix of class "dgCMatrix"
#
#[1,] 1 1 1 . . . . .
#[2,] 1 1 1 . 4 4 . .
#[3,] 1 1 1 . 4 4 . .
#[4,] . . . . 4 4 . .
#[5,] . . 3 3 . . 7 7
#[6,] 2 . 3 3 . . 7 7
#[7,] 2 . . . 5 . . .
#[8,] 2 . . . . 6 6 6
此外,两个函数 return 以相同的顺序分组很方便,我们可以检查:
identical(mySolution(m), ff(m))
#[1] TRUE
举一个看似更复杂的例子:
mm = new("ngCMatrix"
, i = c(25L, 26L, 27L, 25L, 29L, 25L, 25L, 17L, 18L, 26L, 3L, 4L, 5L,
14L, 17L, 18L, 25L, 27L, 3L, 4L, 5L, 17L, 18L, 23L, 26L, 3L,
4L, 5L, 10L, 17L, 18L, 9L, 11L, 17L, 18L, 10L, 17L, 18L, 3L,
17L, 18L, 21L, 17L, 18L, 17L, 18L, 1L, 2L, 3L, 4L, 16L, 8L, 17L,
18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 7L, 9L, 10L, 11L, 26L,
8L, 27L, 1L, 2L, 28L, 1L, 2L, 15L, 27L, 1L, 2L, 21L, 22L, 1L,
2L, 7L, 21L, 22L, 1L, 2L, 6L, 24L, 1L, 2L, 5L, 11L, 16L, 25L,
26L, 27L, 4L, 15L, 17L, 19L, 25L, 26L, 27L, 3L, 16L, 25L, 26L,
27L, 2L, 28L, 1L)
, p = c(0L, 0L, 3L, 3L, 5L, 6L, 7L, 7L, 10L, 18L, 25L, 31L, 35L, 38L,
42L, 44L, 46L, 51L, 61L, 66L, 68L, 71L, 75L, 79L, 84L, 88L, 96L,
103L, 108L, 110L, 111L)
, Dim = c(30L, 30L)
, Dimnames = list(NULL, NULL)
, factors = list()
)
identical(mySolution(mm), ff(mm))
#[1] TRUE
以及在更大矩阵上的简单基准:
times = 30 # times `dim(mm)`
MM2 = do.call(cbind, rep_len(list(do.call(rbind, rep_len(list(mm), times))), times))
dim(MM2)
#[1] 900 900
system.time({ ans1 = mySolution(MM2) })
# user system elapsed
# 449.50 0.53 463.26
system.time({ ans2 = ff(MM2) })
# user system elapsed
# 0.51 0.00 0.52
identical(ans1, ans2)
#[1] TRUE
我的问题在概念上很简单。 我正在寻找一种计算效率高的解决方案(我自己的解决方案附在最后)。
假设我们有一个潜在的非常大的稀疏矩阵,就像下面左边的那个,并且想要 'name' 每个区域的连续非零元素都有一个单独的代码(见右边的矩阵)
1 1 1 . . . . . 1 1 1 . . . . .
1 1 1 . 1 1 . . 1 1 1 . 4 4 . .
1 1 1 . 1 1 . . 1 1 1 . 4 4 . .
. . . . 1 1 . . ---> . . . . 4 4 . .
. . 1 1 . . 1 1 . . 3 3 . . 7 7
1 . 1 1 . . 1 1 2 . 3 3 . . 7 7
1 . . . 1 . . . 2 . . . 5 . . .
1 . . . . 1 1 1 2 . . . . 6 6 6
在我的应用程序中,连续元素将形成矩形、直线或单个点,它们只能通过顶点相互接触(即矩阵中没有 irregular/non 矩形区域)。
我想象的解决方案是将稀疏矩阵表示的行和列索引与具有适当值的向量('name'代码)匹配。我的解决方案使用多个 for loops
并且适用于中小型矩阵,但随着矩阵的维度变大(> 1000),很快就会陷入循环。这可能取决于我在 R 编程方面不是那么先进的事实 - 我找不到任何计算 trick/function 来更好地解决它。
任何人都可以建议一种在 R 中计算更有效的方法吗?
我的解决方案:
mySolution <- function(X){
if (class(X) != "ngCMatrix") {stop("Input must be a Sparse Matrix")}
ind <- which(X == TRUE, arr.ind = TRUE)
r <- ind[,1]
c <- ind[,2]
lr <- nrow(ind)
for (i in 1:lr) {
if(i == 1) {bk <- 1}
else {
if (r[i]-r[i-1] == 1){bk <- c(bk, bk[i-1])}
else {bk <- c(bk, bk[i-1]+1)}
}
}
for (LOOP in 1:(lr-1)) {
tr <- r[LOOP]
tc <- c[LOOP]
for (j in (LOOP+1):lr){
if (r[j] == tr) {
if(c[j] == tc + 1) {bk[j] <- bk[LOOP]}
}
}
}
val <- unique(bk)
for (k in 1:lr){
bk[k] <- which(val==bk[k])
}
return(sparseMatrix(i = r, j = c, x = bk))
}
在此先感谢您的帮助或指点。
在很大程度上依赖于所有要分组的相邻元素仅形成 rectangles/lines/points 这一事实,我们看到矩阵的元素可以根据它们在矩阵上的 [row, col]
索引通过关系聚合(abs(row1 - row2) + abs(col1 - col2)) < 2
。
因此,从 [row, col]
索引开始:
sm = as.matrix(summary(m))
我们计算它们的距离,正如 GiuGe 所说-实际上是 "manhattan" 方法:
d = dist(sm, "manhattan")
Single-linkage 的 属性 在其最近邻元素上的聚类元素在这里很有用。此外,我们可以通过 cutree
ing on "h = 1"(其中索引的距离为“<2”)获得一组元素:
gr = cutree(hclust(d, "single"), h = 1)
最后,我们可以将上面的内容包装在一个新的稀疏矩阵中:
sparseMatrix(i = sm[, "i"], j = sm[, "j"], x = gr)
#8 x 8 sparse Matrix of class "dgCMatrix"
#
#[1,] 1 1 1 . . . . .
#[2,] 1 1 1 . 4 4 . .
#[3,] 1 1 1 . 4 4 . .
#[4,] . . . . 4 4 . .
#[5,] . . 3 3 . . 7 7
#[6,] 2 . 3 3 . . 7 7
#[7,] 2 . . . 5 . . .
#[8,] 2 . . . . 6 6 6
"m"使用的是:
library(Matrix)
m = new("ngCMatrix"
, i = c(0L, 1L, 2L, 5L, 6L, 7L, 0L, 1L, 2L, 0L, 1L, 2L, 4L, 5L, 4L,
5L, 1L, 2L, 3L, 6L, 1L, 2L, 3L, 7L, 4L, 5L, 7L, 4L, 5L, 7L)
, p = c(0L, 6L, 9L, 14L, 16L, 20L, 24L, 27L, 30L)
, Dim = c(8L, 8L)
, Dimnames = list(NULL, NULL)
, factors = list()
)
编辑 2017 年 2 月 10 日
另一个想法(并且再次考虑到相邻元素仅形成 rectangles/lines/points 的事实)是迭代 - 在升序列中 - 通过 [row, col]
索引,并在每一步中找到当前列和行中其最近邻居的每个元素的距离。如果找到 "< 2" 距离,则该元素与其相邻元素分组,否则开始一个新组。包装成一个函数:
ff = function(x)
{
sm = as.matrix(summary(x))
gr = integer(nrow(sm)); ngr = 0L ; gr[1] = ngr
lastSeenRow = integer(nrow(x))
lastSeenCol = integer(ncol(x))
for(k in 1:nrow(sm)) {
kr = sm[k, 1]; kc = sm[k, 2]
i = lastSeenRow[kr]
j = lastSeenCol[kc]
if(i && (abs(kc - sm[i, 2]) == 1)) gr[k] = gr[i]
else if(j && (abs(kr - sm[j, 1]) == 1)) gr[k] = gr[j]
else { ngr = ngr + 1L; gr[k] = ngr }
lastSeenRow[kr] = k
lastSeenCol[kc] = k
}
sparseMatrix(i = sm[, "i"], j = sm[, "j"], x = gr)
}
并应用于 "m":
ff(m)
#8 x 8 sparse Matrix of class "dgCMatrix"
#
#[1,] 1 1 1 . . . . .
#[2,] 1 1 1 . 4 4 . .
#[3,] 1 1 1 . 4 4 . .
#[4,] . . . . 4 4 . .
#[5,] . . 3 3 . . 7 7
#[6,] 2 . 3 3 . . 7 7
#[7,] 2 . . . 5 . . .
#[8,] 2 . . . . 6 6 6
此外,两个函数 return 以相同的顺序分组很方便,我们可以检查:
identical(mySolution(m), ff(m))
#[1] TRUE
举一个看似更复杂的例子:
mm = new("ngCMatrix"
, i = c(25L, 26L, 27L, 25L, 29L, 25L, 25L, 17L, 18L, 26L, 3L, 4L, 5L,
14L, 17L, 18L, 25L, 27L, 3L, 4L, 5L, 17L, 18L, 23L, 26L, 3L,
4L, 5L, 10L, 17L, 18L, 9L, 11L, 17L, 18L, 10L, 17L, 18L, 3L,
17L, 18L, 21L, 17L, 18L, 17L, 18L, 1L, 2L, 3L, 4L, 16L, 8L, 17L,
18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 7L, 9L, 10L, 11L, 26L,
8L, 27L, 1L, 2L, 28L, 1L, 2L, 15L, 27L, 1L, 2L, 21L, 22L, 1L,
2L, 7L, 21L, 22L, 1L, 2L, 6L, 24L, 1L, 2L, 5L, 11L, 16L, 25L,
26L, 27L, 4L, 15L, 17L, 19L, 25L, 26L, 27L, 3L, 16L, 25L, 26L,
27L, 2L, 28L, 1L)
, p = c(0L, 0L, 3L, 3L, 5L, 6L, 7L, 7L, 10L, 18L, 25L, 31L, 35L, 38L,
42L, 44L, 46L, 51L, 61L, 66L, 68L, 71L, 75L, 79L, 84L, 88L, 96L,
103L, 108L, 110L, 111L)
, Dim = c(30L, 30L)
, Dimnames = list(NULL, NULL)
, factors = list()
)
identical(mySolution(mm), ff(mm))
#[1] TRUE
以及在更大矩阵上的简单基准:
times = 30 # times `dim(mm)`
MM2 = do.call(cbind, rep_len(list(do.call(rbind, rep_len(list(mm), times))), times))
dim(MM2)
#[1] 900 900
system.time({ ans1 = mySolution(MM2) })
# user system elapsed
# 449.50 0.53 463.26
system.time({ ans2 = ff(MM2) })
# user system elapsed
# 0.51 0.00 0.52
identical(ans1, ans2)
#[1] TRUE