For prime graph decomposition, one repeatedly looks for cliques that separate the graph.
Please refer to Appendix B.2.2 (page 124) of https://web.math.ku.dk/~lauritzen/papers/gmnotes.pdf for a formal definition of prime graph decomposition.
I am wondering if it has already been implemented in igraph? If not, can I contribute to it? I have already written a function in R that takes an igraph object and output a list of prime components and a list of separators.
Thanks!
Contributions are welcome!
I’m sympathetic towards the idea of including this in igraph, but keep in mind that the igraph team is quite busy, so it’ll take some time to get there.
If we incorporate it into igraph, it should be implemented in C. An R prototype is a good first step. Would you mind sharing your R implementation? Are you comfortable programming in C?
Let’s use this thread to discuss the details.
Thank you! At present I have the following code. I just run through all cliques until I found one that separates the graph. I do not have much experience in C but I don’t mind it.
prime_decom <- function(graph) {
prime_components <- list()
min_sep <- list()
subgraph_stack <- groups(components(graph))
while (length(subgraph_stack) >= 1) {
found = FALSE # a flag
# pop the last graph from stack
curr_graph_v <- subgraph_stack[[length(subgraph_stack)]] # get
subgraph_stack <- subgraph_stack[-length(subgraph_stack)] # remove
curr_graph_v <- sort(curr_graph_v)
curr_graph <- induced_subgraph(graph, curr_graph_v, impl = c("auto"))
all_cliques <- cliques(curr_graph)
for (cli in all_cliques) {
# obtain the induced subgraph by removing cli
remaining_v <- setdiff(V(curr_graph), cli)
removed <- induced_subgraph(curr_graph, remaining_v, impl = c("auto"))
# get connected components of the induced subgraph
clu <- components(removed)
if (clu$no > 1) {
found = TRUE
connected_components <- groups(clu)
# append subgraph_stack
length(subgraph_stack) <- length(subgraph_stack) + clu$no
for (i in 1:clu$no) {
vertices <- union(remaining_v[connected_components[[i]]], cli)
subgraph_stack[[length(subgraph_stack) - i + 1]] <- curr_graph_v[vertices]
}
# append min_sep
length(min_sep) <- length(min_sep) + (clu$no - 1)
for (i in 1:(clu$no - 1)) {
min_sep[[length(min_sep) - i + 1]] <- curr_graph_v[cli]
}
# break the for loop
break
}
}
if (!found) {
# it is a prime component
prime_components <- c(prime_components, list(curr_graph_v))
}
}
return (list(components = prime_components, separators = min_sep))
}
Are you aware of any work on efficient algorithms for decomposition?
Would it work (better) to look for minimal separators and check if they’re cliques?
Have a look at:
https://igraph.org/c/html/latest/igraph-Separators.html#igraph_all_minimal_st_separators
https://igraph.org/c/html/latest/igraph-Separators.html#igraph_minimum_size_separators
Looking for minimal separators is indeed faster. Here I have updated the code.
Since I cannot find the function is_clique in R, I created a subgraph and check the number of edges. I see that this function is available in C.
library(igraph)
prime_decom2 <- function(graph) {
prime_components <- list()
min_sep <- list()
subgraph_stack <- groups(components(graph))
while (length(subgraph_stack) >= 1) {
found = FALSE # a flag
# pop the last graph from stack
curr_graph_v <- subgraph_stack[[length(subgraph_stack)]] # get
subgraph_stack <- subgraph_stack[-length(subgraph_stack)] # remove
curr_graph <- induced_subgraph(graph, curr_graph_v, impl = c("auto"))
## if curr_graph is complete, then prime component
num_edges <- length(curr_graph_v) * (length(curr_graph_v) - 1) / 2
if (ecount(curr_graph) != num_edges) {
all_separators <- min_separators(curr_graph)
for (sep in all_separators) {
induced_sep <- induced_subgraph(curr_graph, sep, impl = c("auto"))
num_edges <- length(sep) * (length(sep) - 1) / 2
if (ecount(induced_sep) == num_edges) {
found = TRUE
curr_graph_v <- sort(curr_graph_v)
# obtain the induced subgraph by removing sep
remaining_v <- setdiff(V(curr_graph), sep)
removed <- induced_subgraph(curr_graph, remaining_v, impl = c("auto"))
# get connected components of the induced subgraph
clu <- components(removed)
connected_components <- groups(clu)
# append subgraph_stack
length(subgraph_stack) <- length(subgraph_stack) + clu$no
for (i in 1:clu$no) {
vertices <- union(remaining_v[connected_components[[i]]], sep)
subgraph_stack[[length(subgraph_stack) - i + 1]] <- curr_graph_v[vertices]
}
# append min_sep
length(min_sep) <- length(min_sep) + (clu$no - 1)
for (i in 1:(clu$no - 1)) {
min_sep[[length(min_sep) - i + 1]] <- curr_graph_v[sep]
}
# break the for loop
break
}
}
}
if (!found) {
# it is a prime component
prime_components <- c(prime_components, list(curr_graph_v))
}
}
return (list(components = prime_components, separators = min_sep))
}
Looks like there is an O(mn) algorithm for this decomposition. I will update my code here once I have it written.
Great find! I would suggest doing a bit of additional literature search before embarking on the implementation. It is useful to find out:
Instead of going for the best known complexity algorithm, at first, I would go for algorithms which seem easier to implement and have a clear description available.
A good next step in literature search is to see what papers cite this one. I haven’t looked at the contents, but the title of this one is promising:
Thank you szhorvat.
I implemented an algorithm presented in a recent paper https://epubs.siam.org/doi/10.1137/16M1059837.
The algorithm goes like this:
(1) Given a graph G, get a minimal chordal completion H.
(2) Get a clique tree of H (so we need the cliques and separators of H, see chordal_factor() below).
[Fact: the minimal separators of H which forms a clique in G are exactly the minimal separators of G.]
(3) Check the minimal separators of H one by one. If the separator does not form a clique in G, then contract that edge in the clique tree.
(4) The resulting clique tree becomes the atom tree, which gives the prime components.
In (1), I rely on the assumption that the output newgraph of the igraph function is_chordal() gives a minimal chordal completion. Could you please verify this?
The algorithm theoretically takes O(n^alpha log n) time, where alpha < 2.3729, using fast matrix multiplication to compute the matrix C (middle of prime_decom3()). Here I simply use %*% in R.
###### input graph must be chordal
chordal_factor <- function(graph) {
n <- vcount(graph)
if (n == 0) {
return(list("cliques" = list(), "separators" = list()))
}
PEO <- as.vector(max_cardinality(graph)$alpham1)
cliques <- list()
kdeg <- vector()
representatives <- vector()
perfectCliques <- list()
separators <- list()
# all possible cliques
for (i in 1:n) {
v <- PEO[i]
neighbours <- intersect(neighbors(graph, v), PEO[i:n])
cliques[[i]] <- c(v, neighbours)
kdeg[i] <- length(cliques[[i]]) - 1
}
# perfect sequence of cliques
perfectCliques <- list(cliques[[1]])
representatives[1] <- PEO[1]
if (n > 1) {
for (i in 2:n) {
if (kdeg[[i]] >= kdeg[[i-1]]) {
representatives <- c(representatives, PEO[i])
perfectCliques <- append(perfectCliques, list(cliques[[i]]))
}
}
perfectCliques <- rev(perfectCliques)
}
# separators
tempUnion <- perfectCliques[[1]]
num_maximal_cliques <- length(perfectCliques)
if (num_maximal_cliques > 1) {
for (i in 2:num_maximal_cliques) {
separators <- append(separators, list(intersect(tempUnion, perfectCliques[[i]])))
tempUnion <- union(tempUnion, perfectCliques[[i]])
}
}
return(list("cliques" = perfectCliques, "separators" = separators))
}
###### input graph is an undirected graph
# assume is_chordal() gives minimal chordal completion
prime_decom3 <- function(graph) {
n <- vcount(graph)
# get a minimal chordal completion of graph
chordal <- is_chordal(graph, newgraph = TRUE)
if (chordal$chordal) {
return (chordal_factor(graph))
} else {
chordal <- chordal$newgraph
factor <- chordal_factor(chordal)
cliques <- factor$cliques
separators <- factor$separators
}
min_sep <- list()
# define clique_matrix and C
clique_matrix <- matrix(0, nrow = n, ncol = length(separators))
for (i in seq_along(separators)) {
clique_matrix[separators[[i]], i] <- 1
}
C <- as_adjacency_matrix(graph, "both") %*% clique_matrix # use fast matrix multiplication
# check
remove <- vector()
for (i in seq_along(separators)) {
sep <- separators[[i]]
if (all(C[as.vector(sep), i] == length(sep) - 1)) {
min_sep <- append(min_sep, list(sep))
} else {
remove <- c(remove, i)
for (j in 1:i) {
if (all(sep %in% cliques[[j]])) {
cliques[[j]] <- union(cliques[[i]], cliques[[j]])
cliques[[i]] <- vector()
break
}
}
}
}
return (list(components = cliques[-remove], separators = min_sep))
}
If is_chordal() does not give minimal chordal completion, I can also implement the algorithm described in Section 4 of the review paper you posted. It’s self-contained.
It turns out that is_chordal() does not give a minimal chordal completion. I found a counter example:
set.seed(124)
graph <- sample_grg(20, 0.4)
chordal <- is_chordal(graph, newgraph = TRUE)$newgraph
print(ecount(chordal)) # 59
print(ecount(graph)) # 49
print(is_chordal(delete_edges(chordal, 54))) # true
I think it is a good idea to modify the implementation of is_chordal() as well (or add a new function for minimal chordal completion if the current chordal completion has other nice properties).
I’m sorry, I forgot to respond here. I suspected that it doesn’t give a minimal one, I didn’t have time to find a reference or example yesterday.
The current algorithm is quite efficient. I’m in favour of writing a new function for a minimal chordal completion.