cluster_leading_eigen() output is incompatible with dendrogram construction

[szhorvat]

Describe the bug

This was reported at https://stackoverflow.com/q/76226275/695132

cluster_leading_eigen() produces a merges matrix in which the indices refer to clusters, not vertices. However, functions that convert to dendrograms still assume that these are vertex IDs, leading to wrong output.

The reason for this sort of output from cluster_leading_eigen() is that the algorithm generates splits (not merges) and produces a dendrogram that is incomplete from below (not incomplete from above, like in the case when we use other algorithms on disconnected graphs).

It seems to me that since the interpretation of the merges matrix returned by cluster_leading_eigen() is different than in the case of other community detection functions, it would be useful not to call it merges (also in the C core) and not to store it in $merges.

To reproduce

g <- make_graph(c(1,5, 1,6, 6,5, 5,4, 2,4, 4,3, 2,3), directed=F)
cl <- cluster_leading_eigen(g)
plot_dendrogram(cl)

Notice the incorrect dendrogram.

For reference, these are the produced clusters:

> groups(cl)
$`1`
[1] 1 5 6

$`2`
[1] 2 3 4

This is the produced merges matrix:

> merges(cl)
     [,1] [,2]
[1,]    1    2

Here 1 and 2 are community indices, not vertex indices.

Additional information

To further illustrate what it means for the dendrogram to be incomplete from below, here's the visualization I use in Mathematica:

It's essentially a dendrogram of communities, not of individual vertices.

Version information

R version 4.2.2 (2022-10-31)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Mojave 10.14.6

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] igraph_1.4.99.9005

loaded via a namespace (and not attached):
 [1] compiler_4.2.2  magrittr_2.0.3  cli_3.6.1       parallel_4.2.2 
 [5] tools_4.2.2     rstudioapi_0.14 Rcpp_1.0.10     nlme_3.1-162   
 [9] ape_5.7-1       grid_4.2.2      digest_0.6.31   rlang_1.1.0    
[13] pkgconfig_2.0.3 lattice_0.21-8 

[szhorvat]

On a related note, this will cause a crash:

cl$merges <-cl$merges + vcount(g)
plot_dendrogram(cl)

[szhorvat]

@ntamas Can you check if python-igraph has a similar issue?

[ntamas]

python-igraph ignores the merge returned from igraph_community_leading_eigenvector() and uses the membership vector and the modularity only, so no, it does not suffer from this problem. But it's not accessible to access the dendrogram either as the Python side does not have a data structure to represent a top-down dendrogram at all.

So yes, I also think that the $merges argument should be renamed to avoid confusion - but it's somewhat hard to do given all the backward compatibilty guarantees that we need to maintain. First of all, what shall we call this member at all? splits?

[szhorvat]

I thought of splits too, as the C documentation suggests it, but it's not really accurate. The difference isn't in splitting vs merging but in the interpretation of the indices.

How about cluster.merges? After all, the indices are cluster indices now.

[szhorvat]

@krlmlr Would you have time to look into the crash I mentioned above? Although to be fair, this is due to corrupting a data structure intentionally. A crash could easily be triggered by messing with the internals of a graph object.

Having ASan set up for CI would make debugging these so much easier ... I usually need 30-60 minutes to set it up locally.

[vtraag]

How about cluster.merges

This sounds like a good suggestion to me.

[szhorvat]

This single change of merges -> cluster.merges would be an easy quickfix that simply disables all hierarchical clustering functionality for cluster_leading_eigen(). Note that is_hierarchical() simply checks for the presence of merges.

At a later stage it would be nice to re-implement this functionality properly for this special case. It's a lot of work though, IMO it's low priority.

In the long term we can think about a unified representation of incomplete dendrograms, and support both those that are "incomplete from below" and "incomplete from above". This might involve the C core as well. Does anyone know if the phylo class from ape supports both?

[ntamas]

I'm not familar with ape but I'm pretty sure that there must be something out in the R ecosystem that supports both types of incomplete dendrograms so we don't need to reinvent the wheel. @krlmlr Any advice here? (i.e. which is the de facto most adopted package in the R ecosystem for dendrograms?).

[krlmlr]

I've only heard about https://talgalili.github.io/dendextend/, not familiar with the packages in that space. @maelle: Can you help?

Regarding the crash, it happens when we call ape::reorder.phylo() in our as_phylo.communities() . Looks like we're constructing a somehow invalid object? I haven't investigated closer.

good_obj <- structure(
  list(
    edge = structure(
      c(7L, 7L, 9L, 9L, 8L, 8L, 11L, 11L, 10L, 10L, 8L, 9L, 10L, 11L, 5L, 6L, 3L, 4L, 1L, 2L),
      dim = c(10L, 2L)
    ),
    edge.length = c(1, 0.5, 1.5, 1, 1.5, 1.5, 1, 1, 0.5, 0.5),
    tip.label = 1:6,
    Nnode = 5L
  ),
  class = "phylo"
)
ape::is.binary(good_obj)
#> [1] TRUE
ape::is.rooted(good_obj)
#> [1] TRUE
ape::reorder.phylo(good_obj)
#> 
#> Phylogenetic tree with 6 tips and 5 internal nodes.
#> 
#> Tip labels:
#>   1, 2, 3, 4, 5, 6
#> 
#> Rooted; includes branch lengths.

bad_obj <- structure(
  list(
    edge = structure(
      c(8L, 8L, 0L, 0L, 10L, 10L, 9L, 9L, 12L, 12L, 11L, 11L, 9L, 10L, 11L, 12L, 5L, 6L, 3L, 4L, 1L, 2L, 7L, 13L),
      dim = c(12L, 2L)
    ),
    edge.length = c(1.5, 1, 2, 1.5, 2, 2, 1.5, 1.5, 1, 1, 0.5, 0),
    tip.label = 1:7,
    Nnode = 6L
  ),
  class = "phylo"
)
ape::is.binary(bad_obj)
#> [1] FALSE
ape::is.rooted(bad_obj)
#> [1] TRUE

# ape::reorder.phylo(bad_obj)

^{Created on 2023-05-20 with reprex v2.0.2}

[maelle]

I don't have any expertise here, I only know of https://docs.ropensci.org/phylogram/