在 R 中获取连通分量

Posted 2023-02-14

【中文标题】在 R 中获取连通分量

【英文标题】：Obtaining connected components in R

【发布时间】：2022-01-06 04:38:42

【问题描述】：

我有一个值为 0 或 1 的矩阵，我想获得一组相邻 1 的列表。

例如矩阵

mat = rbind(c(1,0,0,0,0),
            c(1,0,0,1,0),
            c(0,0,1,0,0),
            c(0,0,0,0,0),
            c(1,1,1,1,1))

> mat
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    0    0    0
[2,]    1    0    0    1    0
[3,]    0    0    1    0    0
[4,]    0    0    0    0    0
[5,]    1    1    1    1    1

应该返回以下 4 个连通分量：

C1 = (1,1);(2,1)

C2 = (2,4)

C3 = (3,3)

C4 = (5,1);(5,2);(5,3);(5,4);(5,5)

有人知道如何在 R 中快速做到这一点吗？我的真实矩阵确实相当大，比如 2000x2000（但我希望连接组件的数量相当少，即 200）。

【参考方案1】：

通过更新，您可以将二进制矩阵转换为栅格对象并使用 clumps 函数。然后它只是返回您想要的确切格式的数据管理。示例如下：

library(igraph)
library(raster)

mat = rbind(c(1,0,0,0,0),
            c(1,0,0,1,0),
            c(0,0,1,0,0),
            c(0,0,0,0,0),
            c(1,1,1,1,1))
Rmat <- raster(mat)
Clumps <- as.matrix(clump(Rmat, directions=4))

#turn the clumps into a list
tot <- max(Clumps, na.rm=TRUE)
res <- vector("list",tot)
for (i in 1:tot)
  res[i] <- list(which(Clumps == i, arr.ind = TRUE))

然后res 在控制台打印出来：

> res
[[1]]
     row col
[1,]   1   1
[2,]   2   1

[[2]]
     row col
[1,]   2   4

[[3]]
     row col
[1,]   3   3

[[4]]
     row col
[1,]   5   1
[2,]   5   2
[3,]   5   3
[4,]   5   4
[5,]   5   5

如果有更好的方法从光栅对象到您的最终目标，我不会感到惊讶。同样，一个 2000 x 2000 矩阵对此应该没什么大不了的。

---

旧的（错误答案），但对于想要连接组件图表的人来说应该很有用。

您可以使用 igraph 包将邻接矩阵转换为网络并返回组件。您的示例图是一个组件，因此我删除了一条边以进行说明。

library(igraph)
mat = rbind(c(1,0,0,0,0),
            c(1,0,0,1,0),
            c(0,0,1,0,0),
            c(0,0,0,0,0),
            c(1,1,1,1,1))
g  <- graph.adjacency(mat) %>% delete_edges("5|3")
plot(g)
clu <- components(g)
groups(clu)

最后一行然后在提示符处返回：

> groups(clu)
$`1`
[1] 1 2 4 5

$`2`
[1] 3

我对这种算法的体验非常快 - 所以我认为 2,000 x 2,000 不会有问题。

【讨论】：

感谢您的回答，但我的矩阵不是邻接矩阵：它不代表图形。我想找到 1 的分组（即矩阵中的 adacent）。

啊抱歉 - connected components 是我所描述的已经使用的术语。看来 raster 包有一个名为 clump 的函数，它可以满足您的需求。我看看能不能举个例子。

太棒了！循环for (i in 1:max(Clumps, na.rm=TRUE)) 可以简化为for (i in 1:tot)，不是吗？