Nonparametric undirected graphical model selection using diffusion models

Technical Summary: Nonparametric Undirected Graphical Model Selection Using Diffusion Models

Problem Statement
Undirected graphical models provide a framework for representing conditional independence structures among high-dimensional random variables. While model selection is a central problem in high-dimensional statistics with applications in genomics, neuroscience, and finance, most existing methods rely on parametric assumptions (e.g., Gaussian graphical models or Ising models). These methods often fail when the true distribution deviates from the assumed parametric form. Nonparametric approaches exist but are statistically challenging due to the "curse of dimensionality" in density estimation and the additional complexity of estimating second-order derivatives (the Hessian of the log-density), which are required to identify the graph structure via the equivalence $(i, j) \notin E_0 \iff \frac{\partial^2 \log p_0}{\partial x_i \partial x_j} = 0$.

Methodology
The authors propose a novel nonparametric approach to undirected graphical model selection based on diffusion models. The method leverages the fact that diffusion models can estimate score functions (gradients of the log-density) and adapt to low-dimensional structures without suffering the full curse of dimensionality.

The proposed procedure, summarized in Algorithm 1, consists of the following steps:

1. Score Function Estimation: Using a standard diffusion modeling framework (specifically the Denoising Diffusion Probabilistic Model, DDPM), the authors estimate a family of score functions $f_0(x, t) = \nabla \log p_t(x)$, where $p_t$ is the marginal density of the forward diffusion process at time $t$. This is achieved via score matching on the observed data.
2. Hessian Estimation via Tweedie's Formula: Instead of directly differentiating the estimated score function (which is numerically unstable), the method utilizes the second-order Tweedie's formula. This formula relates the Hessian of the log-density $\nabla^2 \log p_t(x)$ to the conditional covariance of the original data $X_0$ given the noisy observation $X_t = x$:
   $$\nabla^2 \log p_t(x) = \sigma_t^{-4} \mu_t^2 \text{Cov}[X_0 | X_t = x] - \sigma_t^{-2} I_D$$
   where $\mu_t$ and $\sigma_t$ are parameters of the forward Ornstein-Uhlenbeck process.
3. Graph Recovery: The method estimates the conditional covariance $\text{Cov}[X_0 | X_t = x]$ using Monte Carlo sampling from the reverse diffusion process. Specifically, for a fixed time $t$ and a random vector $U_t$, the authors compute the empirical expectation of the absolute values of the estimated Hessian entries, $E_n[|\hat{H}_{ij}(U_t, t)|]$.
4. Thresholding and Clustering: An edge $(i, j)$ is included in the estimated graph if the averaged magnitude of the Hessian entry exceeds a threshold $\tau$. To handle practical tuning, the authors propose aggregating information across multiple small time steps $t_1, \dots, t_M$ and applying K-means clustering to the standardized Hessian values to separate edges from non-edges.

Key Theoretical Contributions
The paper establishes rigorous theoretical guarantees for the proposed method under mild regularity conditions (specifically, that $\log p_0$ belongs to a Hölder class $\mathcal{H}_{\beta, K}$):

• Decay Rates: The authors prove that the conditional covariance $\text{Cov}[X_{0,i}, X_{0,j} | X_t = x]$ decays at different rates depending on the graph distance $d(i, j)$ between nodes. For edges ($d(i, j)=1$), the covariance decays as $O(t^2)$, whereas for non-edges ($d(i, j) \ge 2$), it decays faster, as $O(t^{(\beta+2)\wedge 5/2})$.
• Selection Consistency: Under the assumption that the score function estimator converges to the true score function at a rate $\epsilon_n$, the authors derive a non-asymptotic upper bound on the probability of incorrect graph selection. They show that with appropriate choices of the time parameter $t$ and threshold $\tau$, the estimator is consistent (i.e., $P(\hat{G}_{t, \tau} \neq G_0) = o(1)$) in a fixed-dimension regime.
• Avoidance of Direct Hessian Estimation: By utilizing the diffusion process and Tweedie's formula, the method circumvents the need to directly estimate the Hessian of the log-density of the original distribution $p_0$, which is statistically more demanding.

Empirical Results
The authors evaluate the method through extensive simulations and two real-data analyses:

• Simulations: The method was tested on four distributions: two Gaussian and two non-Gaussian (Butterfly distribution and Gaussian copula). The proposed diffusion-based method (DDPM) performed competitively with correctly specified parametric models (e.g., Graphical Lasso for Gaussian data) and significantly outperformed existing nonparametric baselines (SING and L-SING) across all distributions. Notably, it achieved perfect graph recovery in non-Gaussian settings where parametric methods failed.
• Image Analysis (MNIST): Applied to the MNIST dataset ($D=784$), the method recovered a sparse conditional independence graph consistent with the locality principle of images (pixels connected primarily to immediate neighbors). It also identified a specific asymmetry in diagonal edges, reflecting the right-handed slant of handwritten digits in the dataset.
• Financial Data (Stock Prices): Applied to daily returns of 28 industrial companies, the method recovered a sparse graph of conditional dependencies. The estimated edges showed meaningful economic relationships, including some not explicitly labeled as competitors in external databases but supported by supply chain or technology partnerships.

Significance and Claims
The paper claims to bridge the gap between the theoretical adaptability of diffusion models to unknown graph structures and their practical application for explicit graph estimation. The authors emphasize that their method offers a unified, nonparametric alternative to existing approaches, capable of handling complex, non-Gaussian distributions without requiring correct model specification. They note that while the theoretical consistency is established for fixed dimensions, the empirical results suggest the method scales effectively to high-dimensional settings ($D \gg n$), leveraging the intrinsic low-dimensional structure of the data. The work concludes by identifying the extension to high-dimensional asymptotics and discrete distributions as natural directions for future research.