Scaling It Up: Stochastic Search Structure Learning in Graphical Models

1 Introduction

• Bayesian graph methods use priors over positive definite matrices with fixed zeros, while classical sparsity methods include penalized likelihood and covariance thresholding.Thresholding does not guarantee that the resulting covariance estimator is positive definite.Source: §1 · pp. 1–2
• Existing Bayesian methods for concentration and covariance graph learning face serious scalability problems, with prior studies reporting long runtimes or only small problems.Reported examples include days for 48- and 100-node problems and 1–20 seconds per one-edge update for a 150-node problem.Source: §1 · p. 2
• SSSL uses continuous shrinkage priors with binary latent indicators to avoid normalizing-constant approximation and enable block graph updates.This contrasts with the point-mass priors used by essentially all existing Bayesian methods for these graph models.Source: §1 · p. 2
• The proposed graph-learning problem differs from regression, factor analysis, and covariance estimation that does not learn graph structure.The paper notes that continuous shrinkage priors had been used in those related problems, while their performance in graph structure learning was little known.Source: §1 · p. 2
• The paper develops two continuous shrinkage prior classes and two block Gibbs samplers for learning undirected and bi-directed graphs.The samplers fit problems with one or two hundred variables within a few minutes.Source: §1 · p. 3

2 Background on graphical models

• Under the multivariate Gaussian model, concentration graphs encode conditional dependence through zeros in Ω, whereas covariance graphs encode marginal dependence through zeros in Σ.The data are modeled as independent samples with covariance matrix Σ, and Ω is its inverse.Source: §2 · p. 3, §2.1 · p. 3, §2.2 · p. 5

2.1 Concentration graph models

• Bayesian concentration graph models use hierarchical priors: a conjugate G-Wishart prior Ω ~ WG(b, D) for Ω given graph Z, and independent Bernoulli priors for each edge-inclusion indicator zij with inclusion probability π.Source: §2.1 · p. 3
• Explicit Bayesian concentration-graph learning is slow because it updates one edge at a time and may require approximation of an intractable normalizing constant.The cited Monte Carlo approximation can require O(Mp^4) computation, and reported desktop runtimes for problems near 100 variables are several days.Source: §2.1 · p. 4
• Continuous-prior estimation of Ω enables fast column-wise block Gibbs updates without normalizing-constant evaluation, but does not formally learn graphs.SSSL is motivated by combining this scalability with explicit graph learning.Source: §2.1 · pp. 4–5

2.2 Covariance graph models

• Covariance graph Bayesian learning uses a G-inverse Wishart prior, but its normalizing constant is analytically intractable and Monte Carlo fitting becomes infeasible beyond a few dozen variables.Closed-form constants for decomposable homogeneous graphs apply to only a small part of the graph space.Source: §2.2 · p. 5
• Covariance graphical lasso and thresholding are alternative approaches, while a block Gibbs fit for a continuous-prior covariance model estimates Σ efficiently without explicit graph learning.The paper proposes combining that computational strength with explicit structure learning.Source: §2.2 · pp. 5–6

3 Continuous spike and slab priors for positive definite matrices

• SSSL places a two-component normal-mixture prior on off-diagonal matrix entries and exponential priors on diagonal entries, normalized over positive-definite matrices.The mixture components use small and large variances, while the normalizing constant ensures the prior integrates to one over the positive-definite matrix space.Source: §3 · p. 6
• Binary latent indicators Z select concentrated or diffuse off-diagonal components, serving as edge indicators whose posterior supports graph inference.An indicator of zero corresponds to a coefficient likely near zero; an indicator of one corresponds to a coefficient that can be substantially nonzero.Source: §3 · p. 7

3.1 Choice of π

• The positive-definiteness constraint biases implied edge inclusion below π, increasingly at larger h, v0, or π, but the authors describe the bias as gentle and data-responsive under a strong likelihood.Figure 1 compares implied inclusion probabilities with reference curves; the reported downward bias does not make the probabilities extremely small.Source: Figure 1 · p. 8, §3.1 · pp. 8–9
• Prior simulation and interpolation can calibrate π to a desired marginal inclusion probability without repeatedly evaluating the intractable normalizer.For p = 150 and target inclusion probability 0.013, the suggested π values are approximately 0.018, 0.027, and 0.048 for v0 = 0.02, 0.05, and 0.1, respectively.Source: §3.1 · p. 9

3.2 Choice of v0 and v1

• Choose v0 so practically negligible values favor zij = 0 and v1 so supported edges can remain substantially nonzero; larger h sharpens the conditional-density separation.The implied conditional densities are not analytically tractable, so threshold-based v0 calibration requires numerical density estimation.Source: §3.2 · pp. 9–10
• For standardized data, MCMC converges quickly and mixes well when v0 ≥ 0.01 and h ≤ 1000.The authors recommend choosing v0 small, but not so small that it causes convergence issues, and choosing v1 to allow reasonable coefficient values.Source: §3.2 · p. 10

3.3 Choice of λ

• For standardized data, λ = 1 assigns substantial probability across plausible diagonal values, and results are insensitive to values such as λ = 5 and 10.The authors report that the data generally contain sufficient information about diagonal elements.Source: §3.3 · p. 10

4 Fast block Gibbs samplers

• SSSL updates all p(p −1)/2 edge-inclusion indicators simultaneously and avoids Monte Carlo approximation of intractable graph constants required by other methods.These two features are identified as the sources of its computational-time reduction.Source: §4 · p. 10
• Posterior graph sampling alternates draws from p(A | Z, Y) and p(Z | A, Y), updating A column-wise and Z all at once.Source: §4 · p. 10

4.1 Block Gibbs samplers for concentration graph models

• Concentration-graph updates are equivalent to p coupled SSVS regressions; using Ω−1_11 instead of S11 in the ω12 conditional shares information across regressions.Source: §4.1 · p. 11

4.2 Block Gibbs samplers for covariance graph models

• For covariance graphs, any column (row) of Σ has a multivariate-normal conditional, enabling direct column-wise block Gibbs updates.Source: §4.2 · p. 12

5.1 Computational efficiency

• At p=100 and 200, 1,000 iterations take about 2 and 9 minutes for covariance graphs, versus 1.2 and 5 minutes for concentration graphs.For p=100, the median inefficiency factor among elements of Ω is 1; 5,000 samples after 2,000 burn-ins usually take under 10 minutes and yield reliable Monte Carlo errors.Source: §5.1 · p. 13, Figure 3 · p. 13

5.2 Structure learning accuracy

• In Scenario 1, larger v0 lowers both TP and FP, larger π raises both, and larger h lowers both through its relation to practical-zero thresholds and edge-inclusion probability.The experiments vary π, h, and v0 across several prior-belief settings.Source: §5.2 · p. 14
• In the p=12 covariance example, SSSL has approximately the G-inverse Wishart’s TP=7 and FP=0, except when hyperparameters favor overly dense or sparse graphs.The adaptive thresholding benchmark has TP=13 and FP=45.Source: §5.2 · p. 15, Table 2 · p. 15
• In the p=150 synthetic Scenario 2, 13% of nonzero partial correlations and 60% of nonzero correlations are within 0.1; weak signals help explain TP drops as v0 increases.Figure 4 shows the empirical frequency histograms for the two datasets.Source: §5.2 · pp. 15, 17, Figure 4 · p. 16
• About 5 minutes versus 2 hours per 1,000 sweeps, SSSL runs faster than Silva (2013), whose sampler updates edges one at a time rather than using SSSL’s block update of Z.Both methods avoid Monte Carlo approximation of the normalizing constant in this comparison.Source: §5.2 · p. 16
• The G-inverse Wishart benchmark yields TP=3 and FP=3 in the p=150 example; the authors attribute its weak graph learning to a prior that concentrates free covariance elements near zero.Their discussion notes that this can make the posterior graph reflect prior influence rather than the data.Source: §5.2 · pp. 16–17
• Further investigation of the G-inverse Wishart framework is left beyond scope; the authors suggest that allowing its parameter b to depend on the graph may help address the restrictive b>0 assumption.For an empty graph, standard inverse-Wishart theory permits b>2−2p, including negative values.Source: §5.2 · p. 17
• With sensible hyperparameters, SSSL is comparable to Bayesian methods; at p=150 it generally improves on adaptive thresholding for covariance graphs and yields sparser concentration graphs than adaptive lasso.The tables report performance under different hyperparameters for p=150 concentration and covariance examples.Source: §5.2 · pp. 17–18, Table 3 · p. 18, Table 4 · p. 18

• The analysis uses daily five-year CDS spread quotes for 79 North American entities and estimates 89 graphs per type from rolling one-year windows.The data span January 2001 through April 2012; the monthly estimation windows are intended to balance observation counts with time variation.Source: §6 · pp. 18–19
• 40–160 of 3,081 possible edges were estimated over time, with both graph types showing similar sparsity and a rise in edge counts from mid-2007.The timing coincides with the start of the subprime-mortgage crisis; the paper says the increase may indicate greater market integration and systemic risk.Source: Figure 5 · p. 19, §6 · p. 19
• By December 2008, the covariance graph had 128 edges and the concentration graph 135, with no isolated nodes in either graph.In December 2004, the respective graphs had 44 and 55 edges, with 30 and 20 isolated nodes; every pair of concentration-graph nodes was connected by a path in 2008.Source: §6 · pp. 19–20, Figure 6 · p. 20
• The authors suggest AIG’s changed ties may reflect the CDS market treating it as unrelated to peer insurers during its more severe credit crisis.They note that such dependence changes may offer a perspective on CDS-market information efficiency distinct from analyses of changes in spread levels.Source: §6 · p. 21
• AIG’s links shifted from three insurance peers in 2004 to GE and AXP in 2008 in both covariance and concentration subgraphs.The links among MET, CB, and HIG remained in the covariance graph after AIG’s links to those firms disappeared.Source: Figure 7 · p. 21, §6 · p. 21

• The proposed framework achieves graph-estimation accuracy comparable to existing Bayesian methods while fitting 100–250-node problems in minutes rather than days or numeric infeasibility.The approach uses absolutely continuous spike-and-slab priors in place of point-mass mixture priors.Source: §7 · p. 22
• Weaker shrinkage near zero may make SSSL parameter estimation suboptimal, although the authors do not find it problematic for structure learning.The paper notes that using distributions with higher densities near zero could refine parameter estimation while retaining data-augmentation-based computational tractability.Source: §7 · p. 22

• The proof considers the conditional distribution of Ω given the edge-inclusion indicators Z.Source: Proof of Proposition 1 · p. 22
• Under partition (10), it considers the conditional distribution of the last column of Ω.Source: Proof of Proposition 1 · p. 22
• It changes variables from (u, ω22) to (u = ω12, v = ω22 − ω′12Ω⁻¹11ω12), with a constant Jacobian independent of those variables.Source: Proof of Proposition 1 · p. 23

• The proof starts from the conditional posterior distribution of Σ given edge-inclusion indicators Z.Source: Proof of Proposition 2 · p. 23
• Under partition (13), it transforms (σ12, σ22) to (u = σ12, v = σ22 − σ′12Σ⁻¹11σ12) and applies block matrix inversion.Source: Proof of Proposition 2 · p. 23
• After removing constants, the proof rewrites the terms as a function of (u, v) and gives their conditionals, including a generalized inverse Gaussian distribution.The passage identifies Kq as a modified Bessel function of the second kind.Source: Proof of Proposition 2 · pp. 23–24