From 2b94f28bf58183a50a2c51dad240833f64f418f5 Mon Sep 17 00:00:00 2001 From: Maarten Marsman Date: Wed, 5 Aug 2026 01:08:37 +0200 Subject: [PATCH] docs: vignette and README pointer for the fast GGM route Adds vignette("fast-ggm"), covering the Gaussian graphical model configuration built for large graphs: precision_graph_prior = "joint" with update_method = "gibbs". The vignette states what the joint specification changes about the prior over graphs before it reports any timing, since the two routes target different models and their inclusion Bayes factors are not interchangeable. It then gives a 50-variable model run both ways (31.4 s default against 3.1 s), a 200-variable model on the fast route (416.7 s for 19,900 candidate edges), a scaling table from 25 to 200 variables, and the configurations the route covers and does not. Timings are precomputed by data-raw/make-fast-ggm.R and shipped in vignettes/fast-ggm.rds, following data-raw/make-prior-sensitivity.R. Every fit ran on its own with four chains on four threads; the hardware, thread count and build are recorded in the payload and printed by the vignette. README gains a short pointer to the route and the vignette. Six references added to vignettes/refs.bib, DOIs verified against Crossref. --- NEWS.md | 10 ++ README.md | 28 +++ data-raw/make-fast-ggm.R | 245 ++++++++++++++++++++++++++ vignettes/fast-ggm.Rmd | 367 +++++++++++++++++++++++++++++++++++++++ vignettes/fast-ggm.rds | Bin 0 -> 1089 bytes vignettes/refs.bib | 67 +++++++ 6 files changed, 717 insertions(+) create mode 100644 data-raw/make-fast-ggm.R create mode 100644 vignettes/fast-ggm.Rmd create mode 100644 vignettes/fast-ggm.rds diff --git a/NEWS.md b/NEWS.md index 35466d79..93bc1097 100644 --- a/NEWS.md +++ b/NEWS.md @@ -43,6 +43,16 @@ released from this line yet. ## Documentation +* New vignette, `vignette("fast-ggm")`, on the Gaussian graphical model route + built for large graphs: `precision_graph_prior = "joint"` with + `update_method = "gibbs"`. It states what the joint specification changes + about the prior over graphs before it reports any timing, gives a fifty- + variable model run both ways and a two-hundred-variable model run on the fast + route, and lists the configurations the route covers. The timings are + precomputed by `data-raw/make-fast-ggm.R` and shipped, as the prior + sensitivity vignette does; the hardware and thread count are stated with them. + The README gained a short pointer to it. + * `prior_sensitivity_check()`'s Details said the refit gate leans on per-chain verdict agreement and the indicator transition ESS. The gate reads the median split-R-hat over the continuous parameters and over the Rao-Blackwellized diff --git a/README.md b/README.md index e0b3d5c3..639ad2b8 100644 --- a/README.md +++ b/README.md @@ -40,6 +40,34 @@ yielding posterior inclusion probabilities for each edge. `bgm()` can additionally model **community structure**, and `bgmCompare()` can test for **group differences** in individual parameters. +## Large Gaussian graphical models + +Continuous data are where graphs get large, and large graphs are where +runtime becomes the constraint. `bgm()` has a faster route for that +case: a conjugate Gibbs sampler for the Gaussian graphical model, +paired with a precision-graph prior that is normalized once rather than +once per graph. Two arguments select it. + +``` r +fit = bgm(y, + variable_type = "continuous", + precision_graph_prior = "joint", + update_method = "gibbs" +) +``` + +On an Apple M5 Pro, four chains on four threads, at the package's +default chain length: a 50-variable graph in 3.1 seconds against 31.4 +for the default route, and a 200-variable graph, 19,900 candidate +edges, in under 7 minutes. + +The joint specification is a different model, not only a faster sampler: +its prior over graphs is the edge prior reweighted, so inclusion Bayes +factors are not interchangeable between the two routes. The vignette +`vignette("fast-ggm")`, also on the [package +website](https://bayesian-graphical-modelling-lab.github.io/bgms/articles/), +sets out that trade-off, the timings, and where the route applies. + ## Installation Install from CRAN: diff --git a/data-raw/make-fast-ggm.R b/data-raw/make-fast-ggm.R new file mode 100644 index 00000000..f302269c --- /dev/null +++ b/data-raw/make-fast-ggm.R @@ -0,0 +1,245 @@ +# ============================================================================== +# Generator for vignettes/fast-ggm.rds +# ============================================================================== +# The fast-GGM vignette reports wall-clock timings. Timings are machine facts, +# and a fit long enough to time is far too long to run inside R CMD check, so +# the runs happen here, once, and the vignette ships and prints the numbers they +# produced. This follows the pattern of data-raw/make-prior-sensitivity.R. +# +# Run from the package root: +# +# Rscript data-raw/make-fast-ggm.R +# +# Every call below MUST stay identical to the call the vignette displays: that +# displayed call is the claim that these numbers came from it. The seeds are +# fixed, and each fit runs on its own, so the elapsed times are not contaminated +# by a neighbouring job. +# +# Machine budget: the runs use 4 chains on 4 threads (chains = 4, cores = 4), +# which is what the vignette states. Nothing here should be run in parallel with +# anything else. +# +# Total runtime is dominated by the 200-variable fit; the whole script is of the +# order of ten minutes on the machine recorded in the payload. +# ============================================================================== + +library(bgms) + +CHAINS = 4L +CORES = 4L + +# ------------------------------------------------------------------------------ +# A sparse Gaussian graphical model. The graph is Erdos-Renyi with an average +# degree of three; the precision matrix is made diagonally dominant so it is +# positive definite by construction, at any number of variables. +# ------------------------------------------------------------------------------ +simulate_sparse_ggm = function(p, avg_degree = 3, weight = 0.25, seed) { + set.seed(seed) + K = matrix(0, p, p) + pairs = which(upper.tri(K), arr.ind = TRUE) + present = stats::runif(nrow(pairs)) < avg_degree / (p - 1) + K[pairs[present, , drop = FALSE]] = + sample(c(-1, 1), sum(present), replace = TRUE) * weight + K = K + t(K) + diag(K) = rowSums(abs(K)) + 1 + K +} + +simulate_case = function(p, n, seed) { + K = simulate_sparse_ggm(p, seed = seed) + list( + K = K, + y = simulate_mrf( + num_states = n, num_variables = p, pairwise = K, + variable_type = "continuous", seed = seed + ), + present = K[upper.tri(K)] != 0 + ) +} + +edge_probabilities = function(fit) { + pip = extract_posterior_inclusion_probabilities(fit) + pip[upper.tri(pip)] +} + +# Recovery against the graph the data came from, at the conventional cut. The +# payload carries these counts and never the inclusion-probability vector +# itself: at 200 variables that vector alone is 20,000 numbers, and the +# vignette quotes summaries of it rather than plotting it. +recovery = function(fit, present) { + pip = edge_probabilities(fit) + list( + selected = sum(pip > 0.5), + detected = sum(pip[present] > 0.5), + false_positive = sum(pip[!present] > 0.5), + max_rhat = max(summary(fit)$indicator[, "Rhat"], na.rm = TRUE) + ) +} + +payload = list() + +# ------------------------------------------------------------------------------ +# Machine and software record. Every timing in the vignette is a fact about this +# machine and this build, and is worthless without it. +# ------------------------------------------------------------------------------ +payload$machine = list( + cpu = tryCatch( + system("sysctl -n machdep.cpu.brand_string", intern = TRUE), + error = function(e) NA_character_ + ), + cores_total = parallel::detectCores(), + os = paste(Sys.info()[["sysname"]], Sys.info()[["release"]]), + r_version = R.version.string, + platform = R.version$platform, + bgms_version = as.character(utils::packageVersion("bgms")), + chains = CHAINS, + threads = CORES, + date = as.character(Sys.Date()) +) + +# ============================================================================== +# 1. The mid-sized model, run both ways +# ============================================================================== +mid = simulate_case(p = 50, n = 1000, seed = 2026) + +payload$mid = list( + p = 50, n = 1000, edges = sum(mid$present), + pairs = length(mid$present) +) + +message("mid-sized (p = 50): default route") +t_default = system.time( + fit_default <- bgm(mid$y, + variable_type = "continuous", + chains = CHAINS, cores = CORES, seed = 2026, + display_progress = "none", verbose = FALSE + ) +) +payload$mid$default = c( + list(elapsed = unname(t_default[["elapsed"]])), + recovery(fit_default, mid$present) +) + +message("mid-sized (p = 50): fast route") +t_fast = system.time( + fit_fast <- bgm(mid$y, + variable_type = "continuous", + precision_graph_prior = "joint", update_method = "gibbs", + chains = CHAINS, cores = CORES, seed = 2026, + display_progress = "none", verbose = FALSE + ) +) +payload$mid$fast = c( + list(elapsed = unname(t_fast[["elapsed"]])), + recovery(fit_fast, mid$present) +) + +payload$mid$speedup = payload$mid$default$elapsed / payload$mid$fast$elapsed +payload$mid$agreement = sum( + (edge_probabilities(fit_default) > 0.5) == + (edge_probabilities(fit_fast) > 0.5) +) + +# The two routes target different models, so the vignette shows what the choice +# does to the prior the analysis is actually run under. Six variables keeps the +# demonstration in seconds: reading the realized prior of a joint fit needs the +# correction table for that model size, and building one is a job of minutes. +message("realized prior inclusion probability (p = 6)") +small = simulate_case(p = 6, n = 300, seed = 1) +fit_small_joint = bgm(small$y, + variable_type = "continuous", + precision_graph_prior = "joint", update_method = "gibbs", + edge_prior = bernoulli_prior(0.5), + chains = 2, cores = 2, iter = 500, warmup = 500, + seed = 1, display_progress = "none", verbose = FALSE +) +fit_small_hier = bgm(small$y, + variable_type = "continuous", + chains = 2, cores = 2, iter = 500, warmup = 500, + seed = 1, display_progress = "none", verbose = FALSE +) +payload$prior = list( + p = 6, + nominal = 0.5, + joint = extract_prior_inclusion_probabilities(fit_small_joint)[1, 2], + hierarchical = extract_prior_inclusion_probabilities(fit_small_hier)[1, 2] +) + +rm(fit_default, fit_fast, fit_small_joint, fit_small_hier) +gc() + +# ============================================================================== +# 2. How the fast route scales, and the high-dimensional fit +# ============================================================================== +# The 200-variable fit is the vignette's high-dimensional example and also the +# last row of the scaling table, so it is run once and reported twice. +# ============================================================================== +grid = c(25, 50, 100, 150, 200) +scaling = data.frame( + p = grid, pairs = grid * (grid - 1) / 2, + edges = NA_integer_, elapsed = NA_real_ +) + +for(i in seq_along(grid)) { + p = grid[i] + message("scaling: p = ", p) + case = simulate_case(p = p, n = 1000, seed = 2026) + scaling$edges[i] = sum(case$present) + tt = system.time( + fit <- bgm(case$y, + variable_type = "continuous", + precision_graph_prior = "joint", update_method = "gibbs", + chains = CHAINS, cores = CORES, seed = 2026, + display_progress = "none", verbose = FALSE + ) + ) + scaling$elapsed[i] = unname(tt[["elapsed"]]) + + if(p == 200) { + payload$high = c( + list( + p = p, n = 1000, pairs = scaling$pairs[i], + edges = scaling$edges[i], elapsed = scaling$elapsed[i] + ), + recovery(fit, case$present) + ) + # How decided the posterior is across all 19,900 pairs. The strongest + # edges are uninformative to tabulate at this size -- there are hundreds + # of them and they all sit at 1.000 -- so what the vignette reports is + # the shape of the whole distribution. + pip = edge_probabilities(fit) + cuts = c(0, 0.05, 0.25, 0.75, 0.95, 1) + payload$high$profile = data.frame( + inclusion_probability = c( + "below 0.05", "0.05 to 0.25", "0.25 to 0.75", + "0.75 to 0.95", "above 0.95" + ), + pairs = as.integer(table(cut(pip, cuts, include.lowest = TRUE))) + ) + # The first six pairs in the fit's own order, as they print. + payload$high$head = utils::head( + summary(fit)$indicator[, c("mean", "sd", "Rhat"), drop = FALSE], 6 + ) + } + rm(fit, case) + gc() +} +payload$scaling = scaling + +# ------------------------------------------------------------------------------ +# The payload ships with the package, so it holds the numbers the vignette +# quotes and nothing else: no chains, no draws, no per-edge vectors. +# ------------------------------------------------------------------------------ +out = file.path("vignettes", "fast-ggm.rds") +saveRDS(payload, out, version = 2) + +size_kb = file.size(out) / 1024 +cat(sprintf("wrote %s (%.1f KB)\n", out, size_kb)) +if(size_kb >= 100) { + stop( + "The vignette payload must stay in double-digit KB; got ", + round(size_kb, 1), " KB. Thin it further before shipping." + ) +} + +str(payload, max.level = 2) diff --git a/vignettes/fast-ggm.Rmd b/vignettes/fast-ggm.Rmd new file mode 100644 index 00000000..3081e51d --- /dev/null +++ b/vignettes/fast-ggm.Rmd @@ -0,0 +1,367 @@ +--- +title: "Speed and Scale in Gaussian Graphical Models" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Speed and Scale in Gaussian Graphical Models} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +bibliography: refs.bib +csl: apa.csl +link-citations: TRUE +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + fig.width = 7, + fig.height = 4.2 +) +library(bgms) +# Every number in this vignette comes from this object. It was produced by one +# run of data-raw/make-fast-ggm.R, whose calls are the calls shown below. +B = readRDS("fast-ggm.rds") +sec = function(x) sprintf("%.1f", x) +# knitr's inline hook renders counts of this size in scientific notation. +num = function(x) formatC(x, format = "d", big.mark = ",") +``` + +# When runtime is the constraint + +Continuous data give a Gaussian graphical model, and continuous data are where +graphs get large. Two hundred variables is twenty thousand candidate edges, and +an analysis that decides them one by one is a long analysis. For that case +`bgm()` has a faster route. Two arguments select it: + +```{r, eval = FALSE} +fit = bgm(y, + variable_type = "continuous", + precision_graph_prior = "joint", + update_method = "gibbs" +) +``` + +`update_method = "gibbs"` replaces the default No-U-Turn sampler with a +conjugate sampler that redraws a whole column of the precision matrix at a time. +`precision_graph_prior = "joint"` normalizes the prior over the precision matrix +and the graph once, globally, instead of once per graph, which puts every term +of the edge move in closed form. Neither is the default, and the second one +changes the model being fitted. That change comes first below; the timings +follow. + +On an Apple M5 Pro, four chains on four threads, at fifty variables the pair is +worth a factor of `r sprintf("%.0f", B$mid$speedup)`. At two hundred variables +(twenty thousand candidate edges) it decides all of them in +`r sprintf("%.0f", B$high$elapsed / 60)` minutes. + + +# The trade-off + +The two routes do not target the same model, and the difference is in the prior +over graphs. + +`bgm()`'s default, `precision_graph_prior = "hierarchical"`, normalizes the +prior over the precision matrix separately within each graph. Integrating the +precision matrix out then returns the stated edge prior exactly: set +`bernoulli_prior(0.5)` and every edge really does carry prior probability +one half. The price is that each edge move has to evaluate a graph-dependent +normalizing constant. + +The joint specification normalizes once for the whole model. That is what makes +the edge move cheap (the constant is the same on both sides of the comparison +and cancels), but it also means the edge prior is no longer the prior over +graphs. What survives is the edge prior reweighted by each graph's normalizing +constant, and that reweighting favors sparser graphs. On six variables, with a +nominal inclusion probability of one half: + +```{r, eval = FALSE} +extract_prior_inclusion_probabilities(fit)[1, 2] +``` + +```{r, echo = FALSE} +data.frame( + route = c("hierarchical (default)", "joint"), + nominal = c(B$prior$nominal, B$prior$nominal), + realized = round(c(B$prior$hierarchical, B$prior$joint), 3) +) +``` + +The default returns the stated half. The joint specification returns +`r round(B$prior$joint, 3)`, and the gap widens as the model grows. Nothing here +is broken (the joint specification is a coherent model, and the sampler targets +its posterior correctly), but it is a *different* model, so its +inclusion Bayes factors are not interchangeable with the default's. An edge that +is undecided on one route can be evidence of absence on the other, without +either fit being wrong. + +`bgm()` says so when it fits: + +``` +Joint precision-graph specification: the realized edge-inclusion prior is the +edge prior reweighted by the per-graph normalizer, not the nominal edge prior +``` + +The choice is a modelling choice, not a tuning knob. Take the fast route when +the size of the problem is the obstacle, and report which specification the +answer came from. An analysis whose report includes a stated prior inclusion +probability needs the nominal edge prior to be the prior over graphs, and that +is the default route. + +One practical consequence. Reading the realized prior back, which +`extract_prior_inclusion_probabilities()`, `extract_inclusion_bf()` and +`verdicts()` all need, requires a table built by simulation for the model size +at hand. It is computed once and cached on disk, but it is not quick: at six +variables it is a few seconds, and by fifty it is longer than the fit it +describes. Posterior inclusion probabilities, from `summary()` or +`extract_posterior_inclusion_probabilities()`, need nothing of the kind. + + +# A mid-sized model, both ways + +`r B$mid$p` variables, `r B$mid$n` observations, drawn from a sparse Gaussian +graphical model with `r B$mid$edges` edges among `r B$mid$pairs` pairs. Both +fits use the package defaults for length (`r formals(bgm)$warmup` warmup and +`r formals(bgm)$iter` retained iterations) and both run four chains on four +threads. They differ in nothing else. + +```{r, eval = FALSE} +default = bgm(y, + variable_type = "continuous", + chains = 4, cores = 4, seed = 2026 +) + +fast = bgm(y, + variable_type = "continuous", + precision_graph_prior = "joint", + update_method = "gibbs", + chains = 4, cores = 4, seed = 2026 +) +``` + +```{r, echo = FALSE} +data.frame( + route = c("hierarchical + nuts (default)", "joint + gibbs"), + seconds = c(round(B$mid$default$elapsed, 1), round(B$mid$fast$elapsed, 1)), + edges_found = c(B$mid$default$selected, B$mid$fast$selected), + true_edges_found = c(B$mid$default$detected, B$mid$fast$detected), + false_positives = c( + B$mid$default$false_positive, B$mid$fast$false_positive + ), + max_rhat = round(c(B$mid$default$max_rhat, B$mid$fast$max_rhat), 3) +) +``` + +`r sec(B$mid$default$elapsed)` seconds against `r sec(B$mid$fast$elapsed)`, a +factor of `r sprintf("%.0f", B$mid$speedup)`, on four chains on four threads. +Both routes recovered `r B$mid$fast$detected` of the `r B$mid$edges` true edges +at the conventional half cut, and both are converged. + +The two edge sets agree on `r B$mid$agreement` of the `r B$mid$pairs` pairs. +That agreement is an observation about this data set, not a guarantee: the +sparser realized prior of the joint specification will pull verdicts toward +absence, and how far it pulls them depends on how much the data have to say. On +data with weaker signal the two routes will part company sooner. The point of +the table is the left-hand column. + + +# Two hundred variables + +`r B$high$p` variables and `r B$high$n` observations: `r num(B$high$pairs)` +candidate edges, `r B$high$edges` of them real. The fast route only. + +```{r, eval = FALSE} +fit = bgm(y, + variable_type = "continuous", + precision_graph_prior = "joint", + update_method = "gibbs", + chains = 4, cores = 4, seed = 2026 +) +``` + +```{r, echo = FALSE} +data.frame( + variables = B$high$p, + candidate_edges = B$high$pairs, + seconds = round(B$high$elapsed, 1), + edges_found = B$high$selected, + true_edges_found = B$high$detected, + false_positives = B$high$false_positive, + max_rhat = round(B$high$max_rhat, 3) +) +``` + +`r sec(B$high$elapsed)` seconds (about +`r sprintf("%.0f", B$high$elapsed / 60)` minutes) for four chains on four +threads at the package's default length, deciding `r num(B$high$pairs)` edges. +It found `r B$high$detected` of the `r B$high$edges` true edges and +`r B$high$false_positive` pairs that are not edges, which on the +`r num(B$high$pairs - B$high$edges)` non-edges is a false-positive rate of +`r sprintf("%.1f%%", 100 * B$high$false_positive / (B$high$pairs - B$high$edges))`. + +Nothing about reading the fit changes at this size: + +```{r, eval = FALSE} +head(summary(fit)$indicator) +``` + +```{r, echo = FALSE} +round(B$high$head, 3) +``` + +The remaining question, across all `r num(B$high$pairs)` pairs, is how decided +the posterior is. The middle of the range holds `r B$high$profile$pairs[3]` +pairs +(`r sprintf("%.1f%%", 100 * B$high$profile$pairs[3] / B$high$pairs)` of the +candidates), so the pairs that need a second look are a short list rather +than a second analysis: + +```{r, echo = FALSE} +knitr::kable(B$high$profile, col.names = c("inclusion probability", "pairs")) +``` + +The default route was not run at this size, so this vignette makes no claim +about what it would have cost. The measured comparison is the one at fifty +variables above. + + +# How it scales + +Same generator, same `r B$high$n` observations, same four chains on four +threads, the fast route throughout. + +```{r, echo = FALSE} +knitr::kable( + data.frame( + variables = B$scaling$p, + candidate_edges = B$scaling$pairs, + seconds = round(B$scaling$elapsed, 1) + ) +) +``` + +```{r, echo = FALSE, fig.alt = "Wall-clock seconds against the number of variables, on a logarithmic vertical axis, rising steeply from a fraction of a second at 25 variables to several minutes at 200."} +op = par(mar = c(4.2, 4.4, 1.2, 1.2), las = 1, bty = "n", cex = 0.9) +plot(B$scaling$p, B$scaling$elapsed, + log = "y", type = "b", pch = 19, lwd = 2, cex = 1.1, + col = "#2f6f9f", + xlab = "Number of variables", + ylab = "Seconds (four chains, four threads)", + panel.first = grid(col = "grey92", lty = 1, lwd = 1) +) +text(B$scaling$p, B$scaling$elapsed, + labels = sprintf("%.1f", B$scaling$elapsed), + pos = c(4, 4, 4, 4, 2), offset = 0.6, cex = 0.8, col = "grey30" +) +par(op) +``` + +What drives that curve is the number of variables, not the number of +observations. Complete data enter the sampler only through the cross-product of +the data matrix, formed once before sampling starts; every draw after that +works on the precision matrix, whose size the variables alone set. More +observations change what the sampler is drawing from, not how much work a draw +costs. More variables change both. + + +# Where the fast route applies + +`update_method = "gibbs"` is the Gaussian graphical model's sampler and nothing +else. Ordinal, Blume-Capel or mixed data are refused outright: + +```{r, error = TRUE} +bgm(Wenchuan[, 1:4], update_method = "gibbs") +``` + +Within continuous data it needs conjugate ingredients: a `normal_prior()` or +`cauchy_prior()` on the interactions, and `exponential_prior()` or +`gamma_prior()` on the precision diagonal. Those are the defaults for +`interaction_prior` and `precision_scale_prior`, so the usual call already +qualifies; `beta_prime_prior()` on the interactions does not. Edge selection is +supported for both slabs, and so is `edge_selection = FALSE`, which estimates +the precision matrix at a fixed graph. + +`precision_graph_prior = "joint"` is separate from the sampler and applies to +continuous data under edge selection, in Gaussian graphical models and in the +continuous block of a mixed model. The two arguments compose but do not depend +on each other: the Gibbs sampler runs under the default hierarchical +specification, and the joint specification runs under `"nuts"`. + +One combination is fast to fit but slow to start. Under the joint +specification a learned inclusion probability (`beta_bernoulli_prior()` or +`sbm_prior()`) needs the same simulated correction table described above, +because the constant that cancels for a fixed inclusion probability stops +cancelling once that probability is itself being updated. The table is built at +the first fit of a model configuration and cached, but at large numbers of +variables that first fit will be waiting on it for far longer than the sampling +takes. For a fast run at scale, keep `bernoulli_prior()`, the default. + + +# Background + +The prior used here specifies the precision matrix entry by entry: an edge that +is absent has its entry fixed at exactly zero, and an edge that is present draws +its entry from a diffuse slab. The construction comes from the Bayesian +variable-selection literature [@george1993jasa_variable] and was brought to +Gaussian graphical models by @wang2015ba_scaling, whose column-wise block Gibbs +sampler is the ancestor of the one `update_method = "gibbs"` runs. Its appeal +was, and is, that it avoids the graph-dependent normalizing constant that makes +structure learning expensive under the G-Wishart prior +[@atay-kayis2005biometrika_monte; @mohammadi2023jasa_accelerating]. That +avoidance is exactly what the joint specification does and what the +hierarchical specification declines to do, which is the choice +`precision_graph_prior` exposes. @vogels2024jasa_bayesian survey the wider field +and compare the available samplers empirically. + +`bgms` adds one thing to the classical construction. Specifying the precision +matrix entry by entry does not guarantee that the result is positive definite, +so the prior has to be truncated to the cone of positive-definite matrices, and +that truncation distorts it more and more as the number of variables grows +[@jewson2024biometrics_graphical]. `bgm()` counters this by tilting the prior by +a power `delta` of the determinant, which pushes mass away from the boundary of +the cone; `delta` defaults to half the logarithm of the number of variables, so +larger models are tilted harder. This is why the numbers above hold up at two +hundred variables rather than collapsing onto the empty graph. + +Scaling this prior is an active problem and `bgms` is not alone on it. +@sulem2025arxiv_bayesian give the closest treatment: +Metropolis-Hastings within block Gibbs on the same spike-and-slab Gaussian +graphical model, with mixing bounds that are dimension-free under suitable +conditions and a global proposal that adds or removes several edges in one +move. They report carrying exact Bayesian inference from roughly one hundred +variables to roughly one thousand. Their work is implemented in the CRAN package +**modelSelection** (version 1.0.7, 2026-05-16), which covers Bayesian model +selection and averaging for regression, generalized linear and additive models, +mixtures and graphical models. **ssgraph** (version 1.16, 2025-08-29) is the +other near neighbor, doing Bayesian structure learning for undirected graphical +models with spike-and-slab priors on continuous, discrete and mixed data. + +Neither package was run against `bgms` on the same data for this vignette, so it +makes no comparative speed claim about either. Every timing above is `bgms` +against `bgms`. Readers whose problem is a single very large Gaussian graphical +model and nothing else should look at all three; what `bgms` adds around this +route is the rest of the package: ordinal, Blume-Capel and mixed models, +group comparison, prior sensitivity, and the determinant tilt that keeps the +prior usable as the number of variables grows. + + +# Reproducing the timings + +Every number above was produced by `data-raw/make-fast-ggm.R`, which is in the +package repository and shows the generator, the seeds and the calls. Each fit +ran on its own, four chains on four threads. + +```{r, echo = FALSE} +data.frame( + field = c("CPU", "cores available", "threads used", "chains", + "operating system", "R", "bgms", "run"), + value = c(B$machine$cpu, B$machine$cores_total, B$machine$threads, + B$machine$chains, B$machine$os, B$machine$r_version, + B$machine$bgms_version, B$machine$date) +) +``` + +Timings are facts about a machine. The ratios will travel further than the +seconds. + + +# References diff --git a/vignettes/fast-ggm.rds b/vignettes/fast-ggm.rds new file mode 100644 index 0000000000000000000000000000000000000000..ca32f769e5230d90c96ee3af22c23f1359a94d5e GIT binary patch literal 1089 zcmV-H1it$piwFP!000001HDyWj1xr|-)?t%ZI1&DIsfI*g+D{u2FV#QG3u!iA|zZq@Xc@!MhOolC@&b`1(X-j5Q$L}dv|BFJaLMgN47u%x->jR zjq!>zO~w9-jiFoeY!8^!B-$>OQn*CVMGZ;E1YXs~4VJxft>H#u5|uhY7TcoFt78ZZ zrb-6zAesN5*O4S+xevGapum2C{(6z=?9f{}vUI9*tu$UoA5P`y0EIGD|Jk*3NssWxZ;X3fR zLw7(Yul{)S{zGRuaNDi{I?(|weLA?fU;iILQ?SEBkF|bTvy!`o6A=(ex<2MGffa=e zjT&S{AUgozus{r01=FU-3Zx2LCil4xgvK=6ECp&9!UnYXTuLFlo4N}~zZlWaZ~SBu zY@Ye6zv~E6^ay>H-9T2@9P$%3vKYW+bp)~B&Zt5HI-vSPP^cROSJi=j&SpBzX>4`q zd$#|U`Hz|>IUiTuF_9~dot@%HWbR${@IKPwJ`OQb9O*JN?UtgcgeK?-d|B~qV#{*2 zEVYeswlQiOXSOsb+f5i)0xYDk@7Yg2St3x55q+GDGm!*~boC&zD(1A$8zl$ZWl%;Y z@}B9gkO-AkOyDvq3z?4)3Q((ifu(Esgq`?uJYmHun< z4^KH=Wb_{sCf$^jbQPJ*E!D{r$Wfw7rXc}2Hf7d?&T@XCC HB?