如何用R的igraph包高效计算图的双星/楔数量？求优化方案

问题背景 #

igraph包提供了高效的三角形计数函数count_triangles()，但无法通过getAnywhere()获取其源码。我推测基于该函数的思路可以实现高效的双星（楔）计数，但目前自己实现的楔计数方法在大图中性能极差，希望得到优化建议。

示例图与现有实现 #

先定义示例无向图：

w3 <- matrix(c(c(0,1,1,1), c(1,0,1,0), c(1,1,0,0), c(1,0,0,0)), nrow = 4, ncol = 4)
w21 <- graph_from_adjacency_matrix(w3, mode = "undirected")

三角形计数的性能对比 #

我实现了两种三角形计数方法：

# 基于邻接矩阵乘法的方法
triangle_cnt <- function(adj_mat){
  return(sum(diag(adj_mat %*% adj_mat %*% adj_mat)))
}

# 基于igraph的方法
triangle_cnt.igraph <- function(adj_mat){
  adg_mat <- graph_from_adjacency_matrix(adj_mat)
  return(sum(count_triangles(adg_mat)) * 2)
}

# 性能测试
microbenchmark::microbenchmark(triangle_cnt(w3), triangle_cnt.igraph(w3))

小图测试输出：

> microbenchmark::microbenchmark(triangle_cnt(w3), triangle_cnt.igraph(w3))
Unit: microseconds
                    expr    min     lq     mean median      uq     max neval
        triangle_cnt(w3)  1.435  1.640  2.06394   2.05  2.1730  14.104   100
 triangle_cnt.igraph(w3) 57.400 58.999 63.25521  59.86 62.3815 254.692   100

大图（1000+节点）测试输出：

> microbenchmark::microbenchmark(triangle_cnt(adj_mat), triangle_cnt.igraph(adj_mat))
Unit: milliseconds
                         expr      min       lq     mean   median       uq      max neval
        triangle_cnt(adj_mat) 651.1529 660.7670 669.5889 663.1095 665.9754 753.3584   100
 triangle_cnt.igraph(adj_mat) 180.1882 181.3038 195.4430 181.9534 182.6366 524.7692   100

现有楔计数方法的性能问题 #

我当前的楔计数方法基于邻接矩阵相乘求和，在大图中速度极慢：

wedge_cnt <- function(adj_mat){
  return(sum(adj_mat %*% adj_mat))
}

优化方案：利用节点度计算楔数量 #

对于无向图，楔的总数等于所有节点的度的组合数之和，公式为：
$$\text{总楔数} = \sum_{i=1}^n \binom{d_i}{2} = \sum_{i=1}^n \frac{d_i \times (d_i - 1)}{2}$$
其中$d_i$是第$i$个节点的度。这个公式的逻辑是：每个楔由一个中心节点连接两个邻居构成，每个度为$d_i$的节点能贡献$\binom{d_i}{2}$个这样的楔，累加所有节点的贡献即可得到总楔数。

实现代码 #

基于igraph的高效实现

wedge_cnt.igraph <- function(graph){
  degrees <- degree(graph)
  return(sum(degrees * (degrees - 1) / 2))
}

# 用示例图测试
wedge_cnt.igraph(w21)

基于邻接矩阵的轻量实现（无需转换为igraph对象）

wedge_cnt.fast <- function(adj_mat){
  degrees <- rowSums(adj_mat)
  return(sum(degrees * (degrees - 1) / 2))
}

# 用示例图测试
wedge_cnt.fast(w3)

性能优势说明 #

这种方法完全避免了O(n³)的矩阵乘法，仅需O(n)的计算（计算节点度）和O(n)的求和，在大图上的性能会远优于矩阵相乘的方法。比如对于1000节点的图，该方法的计算时间会从数百毫秒级降至微秒级。

内容的提问来源于stack exchange，提问作者ann