Weighted Crossing Minimization: wompwomp vs. Classical Heuristics

Background

The two-layer (“bipartite”) crossing-minimization problem underlies many diagram types, from alluvial/Sankey plots to storyline charts, and is the subject of dedicated benchmarks such as the 2024 PACE Challenge. Two classic, cheap heuristics for that problem are the barycenter heuristic (place each node at the average position of its neighbors) and the median heuristic (place each node at the median position instead) – both available here as method = "barycenter"/"median" in [sort_to_uncross()].

That classical problem is unweighted: every crossing counts the same. Alluvial data is not – edges have a magnitude (the size of the flow), and sort_to_uncross()‘s default objective (weighted_metric = TRUE) weights each crossing by the product of the two crossing edges’ magnitudes, so a crossing between two large flows counts for more than one between two small ones. neighbornet (this package’s default heuristic, adapted from phylogenetic network reconstruction) and the two-sided greedy heuristic were both designed with this weighted objective in mind. This vignette checks whether that design choice actually pays off: how close does each method get to the true weighted-optimal ordering, on data built for the weighted case?

Establishing ground truth

We compare five methods – tsp, neighbornet, greedy, barycenter, and median – against the true optimum. For small enough instances, that optimum can be found by brute force: try every ordering of both layers and keep the best.

combinat_perms <- function(x) {
    if (length(x) == 1) return(list(x))
    res <- list()
    for (i in seq_along(x)) for (r in combinat_perms(x[-i])) res[[length(res) + 1]] <- c(x[i], r)
    res
}

# Exhaustive search over both layers' orderings -- (n_categories!)^2 candidate
# arrangements, so this only stays fast enough to keep here for small
# n_categories (up to ~5).
brute_force_optimum <- function(df, n_categories) {
    labels <- as.character(seq_len(n_categories) - 1L)
    perms <- combinat_perms(labels)
    best <- Inf
    for (pa in perms) {
        for (pb in perms) {
            d2 <- df
            d2$A <- factor(d2$A, levels = pa)
            d2$B <- factor(d2$B, levels = pb)
            obj <- compute_crossing_objective(d2, cols = c("A", "B"), wt = "value", weighted_metric = TRUE)$output_objective
            if (obj < best) best <- obj
        }
    }
    best
}

Synthetic weighted data

Uniform-random categories give every heuristic no structure to exploit, which isn’t representative of real alluvial data (e.g. tissue/cell-type combinations, where columns correlate through a shared underlying identity). This generator gives each row a latent cluster label shared across both columns; with probability cluster_strength a column’s value follows that cluster (through its own random permutation of labels), otherwise it’s uniform noise.

make_weighted_two_layer_df <- function(n_rows, n_categories, seed, cluster_strength = 0.7) {
    set.seed(seed)
    labels <- as.character(seq_len(n_categories) - 1L)
    cluster <- sample.int(n_categories, size = n_rows, replace = TRUE)
    col_perms <- replicate(2, sample(labels), simplify = FALSE)
    raw <- as.data.frame(lapply(1:2, function(j) {
        matches_cluster <- stats::runif(n_rows) < cluster_strength
        noise <- sample(labels, size = n_rows, replace = TRUE)
        ifelse(matches_cluster, col_perms[[j]][cluster], noise)
    }), stringsAsFactors = FALSE)
    names(raw) <- c("A", "B")
    as.data.frame(dplyr::count(raw, A, B, name = "value"))
}

Running the comparison

For a few problem sizes and several random replicates each, find the true optimum and see how far each method’s result lands from it.

methods <- c("tsp", "neighbornet", "greedy", "barycenter", "median")
n_categories_vals <- c(3, 4, 5)
n_seeds <- 6

rows <- list()
for (nc in n_categories_vals) {
    for (seed in seq_len(n_seeds)) {
        df <- make_weighted_two_layer_df(n_rows = 500, n_categories = nc, seed = seed)
        optimum <- brute_force_optimum(df, nc)
        for (method in methods) {
            sorted <- sort_to_uncross(df, cols = c("A", "B"), wt = "value", method = method)
            obj <- compute_crossing_objective(sorted, cols = c("A", "B"), wt = "value", weighted_metric = TRUE)$output_objective
            rows[[length(rows) + 1]] <- data.frame(
                n_categories = nc, seed = seed, method = method,
                objective = obj, optimum = optimum
            )
        }
    }
}
results <- do.call(rbind, rows)
results$pct_above_optimal <- 100 * (results$objective - results$optimum) / results$optimum

Results

summary_table <- aggregate(
    pct_above_optimal ~ method, data = results,
    FUN = function(x) c(median = round(median(x), 2), mean = round(mean(x), 2), max = round(max(x), 2))
)
summary_table <- do.call(data.frame, summary_table)
names(summary_table) <- c("method", "median_pct_above_optimal", "mean_pct_above_optimal", "max_pct_above_optimal")
knitr::kable(summary_table[order(summary_table$median_pct_above_optimal), ])

if (requireNamespace("ggplot2", quietly = TRUE)) {
    results$method <- factor(results$method, levels = summary_table$method[order(summary_table$median_pct_above_optimal)])
    ggplot2::ggplot(results, ggplot2::aes(x = method, y = pct_above_optimal)) +
        ggplot2::geom_boxplot(outlier.shape = NA) +
        ggplot2::geom_jitter(width = 0.15, alpha = 0.5) +
        ggplot2::labs(
            x = NULL, y = "% above true optimal (weighted crossings)",
            title = "Distance from the true weighted-optimal ordering",
            subtitle = sprintf("%d configs (n_categories in %s, %d seeds each)", nrow(results) / length(methods), paste(n_categories_vals, collapse = ", "), n_seeds)
        ) +
        ggplot2::theme_minimal()
}

Discussion

neighbornet – this package’s default two-sided method – lands closest to the true optimum (median 0% above optimal across these 18 configurations, hitting the exact optimum in more than half of them), clearly ahead of the classical barycenter/median heuristics (median ~8.6% above optimal) and the two-sided greedy heuristic, which performs identically to barycenter/median at this scale. That three otherwise-different algorithms converge to the same answer here isn’t a bug: at small n_categories with strongly clustered data, there’s often one dominant, “obvious” arrangement that several reasonable heuristics independently find – confirmed above by checking that they diverge normally once the clustering signal is weakened. tsp performs comparably to neighbornet here, though with a higher ceiling on how far it can miss (up to 15.9% above optimal in this sample, versus 6.8% for neighbornet) despite being exponentially more expensive to run.

This is a small, brute-force-validated demonstration rather than a large-scale benchmark – see benchmarks/ in the package source for larger sweeps – but it directly answers the motivating question: on data built for the weighted, correlated case sort_to_uncross() targets, its default heuristic does outperform the classical unweighted-problem heuristics, not just match them.

sessionInfo()

## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 26.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] wompwomp_1.0.0    ggalluvial_0.12.6 ggplot2_4.0.3     dplyr_1.2.1      
## [5] rmarkdown_2.32   
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.10        utf8_1.2.6         generics_0.1.4     tidyr_1.3.2       
##  [5] stringi_1.8.9      digest_0.6.39      magrittr_2.0.5     evaluate_1.0.5    
##  [9] grid_4.6.1         RColorBrewer_1.1-3 iterators_1.0.14   fastmap_1.2.0     
## [13] foreach_1.5.2      jsonlite_2.0.0     purrr_1.2.2        scales_1.4.0      
## [17] tweenr_2.0.3       codetools_0.2-20   jquerylib_0.1.4    cli_3.6.6         
## [21] rlang_1.3.0        polyclip_1.10-7    withr_3.0.3        cachem_1.1.0      
## [25] yaml_2.3.12        otel_0.2.0         tools_4.6.1        buildtools_1.0.0  
## [29] vctrs_0.7.3        TSP_1.2.7          R6_2.6.1           lifecycle_1.0.5   
## [33] stringr_1.6.0      MASS_7.3-66        pkgconfig_2.0.3    pillar_1.11.1     
## [37] bslib_0.12.0       gtable_0.3.6       glue_1.8.1         Rcpp_1.1.2        
## [41] ggforce_0.5.0      xfun_0.62          tibble_3.3.1       tidyselect_1.2.1  
## [45] sys_3.4.3          knitr_1.52         farver_2.1.2       htmltools_0.5.9   
## [49] igraph_2.3.4       maketools_1.3.2    labeling_0.4.3     compiler_4.6.1    
## [53] S7_0.2.2