Likelihood-Free Inference via Structured Score Matching

1. Problem Statement

The paper addresses the challenge of Likelihood-Free Inference (LFI), also known as Simulation-Based Inference (SBI). In many scientific fields (e.g., astrophysics, epidemiology, genetics), complex phenomena are modeled using forward simulators. While these simulators can generate data $X \sim P_\theta$ for any parameter $\theta$, the likelihood function $p_\theta(X)$ is often analytically intractable or unavailable in closed form.

• Goal: Perform frequentist inference (parameter estimation and uncertainty quantification) without evaluating the likelihood.
• Objective: Estimate the Maximum Likelihood Estimator (MLE), $\hat{\theta}_{MLE}$, and construct valid confidence sets, despite the inability to compute $\nabla_\theta \log p_\theta(X)$.
• Limitations of Existing Methods:
  • Approximate Bayesian Computation (ABC): Relies on summary statistics and distance metrics, often losing statistical efficiency.
  • Neural Likelihood Estimators (NLE): Directly approximate the posterior or likelihood using deep generative models. These can be biased, poorly calibrated in low-density regions, and heavily dependent on training objectives.
  • Indirect Inference: Often relies on auxiliary models that may not fully capture the true data structure.

2. Methodology

The authors propose a framework that approximates the likelihood score function (the gradient of the log-likelihood) rather than the likelihood itself. The core innovation is Structured Score Matching, which embeds fundamental statistical properties of the true score directly into the neural network architecture and training objective.

A. Structured Score Matching

Instead of learning a generic function $s_\phi(\theta, X_n)$ to approximate the full-data score $s^*(\theta, X_n) = \nabla_\theta \sum \log p_\theta(X_i)$, the method enforces three key statistical properties:

1. Additive Structure: The score of independent samples is the sum of individual scores.
   • Implementation: The network is constrained to output $s_\phi(\theta, X_n) = \sum_{i=1}^n s_\phi(\theta, X_i)$. This reduces the input dimension from $n \times p$ to $p$, significantly lowering the complexity and sample complexity of training.
2. Curvature Structure (Fisher Information Identity): The expected outer product of the score plus the expected Jacobian of the score equals zero: $E[s s^\top + \nabla s] = 0$.
   • Implementation: A regularization penalty is added to the loss function to enforce this relationship, ensuring the estimated score respects the local geometry of the log-likelihood.
3. Mean-Zero Structure: The expected score is zero: $E[s] = 0$.
   • Implementation: A two-step debiasing procedure is used. First, a score network is trained. Second, a regression model is trained to estimate the mean bias of the score, which is then subtracted to produce a debiased estimator.

The training objective minimizes the expected squared distance between the estimated score and the true score, utilizing the Hyv¨arinen (2007) identity to avoid computing the intractable true score.

B. Score-Based Estimator

Once the structured score function $\hat{s}(\theta, X_n)$ is trained, the parameter estimator $\hat{\theta}_n$ is defined as the root of the estimated score:
$$\hat{\theta}_n = \arg \min_{\theta} \| \sum_{i=1}^n \hat{s}(\theta, X_i) \|^2$$
To find this root, the authors propose an iterative algorithm (Algorithm 4) that mimics gradient ascent on the log-likelihood:
$$\theta^{(t+1)} = \theta^{(t)} + \alpha H_t \hat{s}(\theta^{(t)}, X_n)$$

• Preconditioning ($H_t$): The method leverages the curvature matching property to use the inverse Jacobian of the estimated score ($H_t \approx [\nabla_\theta \hat{s}]^{-1}$), effectively implementing a Newton-Raphson or Quasi-Newton update. This ensures exponential convergence rates locally, unlike standard gradient descent.

C. Uncertainty Quantification

The paper provides three methods to construct confidence sets for $\theta^*$:

1. Information Matrix Plug-in: Uses the empirical average of the negative Jacobian of the estimated score to approximate the Fisher Information Matrix.
2. Sandwich Covariance (Huber): Treats the estimator as an M-estimator under potential model misspecification. It uses the sandwich formula: $\hat{\Sigma} = \hat{I}^{-1} \hat{J} \hat{I}^{-1}$, where $\hat{J}$ is the empirical variance of the scores. This is shown to be more robust.
3. Multiplier Bootstrap: A resampling method where weighted scores are solved to generate bootstrap distributions, providing asymptotically valid confidence intervals without relying on normality assumptions.

3. Key Contributions

1. Theoretical Guarantees:
   • Proves the existence and uniqueness of the root of the structured score function in a neighborhood of the true parameter.
   • Establishes that the proposed estimator is asymptotically equivalent to the MLE (i.e., $\sqrt{n}(\hat{\theta}_n - \hat{\theta}_{MLE}) \to 0$).
   • Proves asymptotic normality of the estimator and bootstrap consistency.
   • Demonstrates exponential convergence of the iterative optimization algorithm due to the curvature matching.
2. Architectural Innovation:
   • Introduces Structured Score Matching, which explicitly encodes additivity, curvature, and mean-zero properties into the neural network. This improves accuracy and scalability compared to "vanilla" score matching.
3. Two-Round Procedure:
   • Proposes a refinement strategy where an initial coarse estimate is used to define a focused sampling distribution for a second round of training, improving score approximation in the region of interest.

4. Experimental Results

The method was evaluated on synthetic and real-world datasets against Synthetic Likelihood (SL), Neural Likelihood Estimation (NLE), and Iterative Local Score Matching (ILSM).

• M/G/1 Queuing Model: The proposed method achieved lower estimation errors than SL and NLE. The Sandwich and Bootstrap confidence intervals achieved near-nominal coverage (92-94%) with tighter widths compared to SL. The Quasi-Newton optimizer converged in ~4 iterations vs. ~81 for gradient descent.
• Stock Volatility Estimation: On a multivariate Brownian motion model, the method yielded the smallest estimation errors for all parameters. Confidence intervals constructed via plug-in or bootstrap methods had near-nominal coverage and significantly smaller widths than NLE and SL.
• g-and-k Model (Real Data): Applied to US Dollar/Canadian Dollar exchange rates. The method successfully modeled heavy tails and skewness. The two-round procedure showed clear improvements in estimation accuracy and coverage over the first round.
• Efficiency: The method requires fewer simulator calls than SL (which relies on MCMC) and comparable calls to NLE. While training time is higher than SL (due to neural networks), it is competitive with NLE and offers superior statistical properties.

5. Significance

This work bridges the gap between classical frequentist statistics and modern deep learning for simulation-based inference.

• Statistical Efficiency: By targeting the score function (which leads to the MLE) rather than the likelihood or posterior directly, the method retains the asymptotic efficiency of the MLE.
• Robustness: The incorporation of statistical structure (curvature/additivity) prevents the "black box" nature of standard neural approximations, ensuring the estimator behaves correctly as sample size grows.
• Practical Utility: The provision of rigorous uncertainty quantification (confidence intervals) via multiple robust methods (Sandwich, Bootstrap) makes the approach suitable for scientific applications where decision-making requires reliable error bounds.
• Scalability: The additive structure allows the method to scale to larger datasets without retraining the network on full data batches, a significant advantage over methods that treat the full dataset as a single input.

In summary, the paper presents a theoretically grounded, highly efficient, and statistically robust framework for likelihood-free inference that outperforms existing state-of-the-art methods in both accuracy and uncertainty calibration.