Pseudo-Likelihood: Methods & Applications

• Pseudo-likelihood is a surrogate likelihood method that replaces full joint distributions with products of local conditional models, bypassing intractable normalizing constants.
• Under regularity and identifiability conditions, pseudo-likelihood estimators achieve consistency and asymptotic normality, offering practical inference in high-dimensional settings.
• Widely applied in graphical models, inverse-covariance estimation, Bayesian structure learning, and distributed computation, pseudo-likelihood provides computational tractability with theoretical assurances.

Searching arXiv for recent and foundational papers on pseudo-likelihood to ground the article. arXiv search query: pseudo-likelihood graphical models consistency pseudo-likelihood inference Pseudo-likelihood is a surrogate likelihood principle in which an intractable joint model is replaced by a product, or sum of logarithms, of tractable local conditional models or related score components. In the classical form introduced for Markov networks, the global partition function is removed by conditioning on local neighborhoods, so inference proceeds through quantities that normalize only over a single variable or a small block rather than over the full state space (Pensar et al., 2014). In contemporary usage, the term covers a broader family of estimators and scoring rules for graphical models, copula models, stochastic block models, sequence models, dynamic games, and simulation-based inference, where the common aim is to preserve inferential content while avoiding the computational bottlenecks of exact likelihood evaluation (Mozeika et al., 2014).

1. Core construction and relation to exact likelihood

For discrete observations $s=(s_1,\dots,s_N)$, a standard maximum pseudo-likelihood objective is

$$PLL_L(\theta) = \frac{1}{L}\sum_{\mu=1}^L \sum_{i=1}^N \log P_\theta\!\bigl(s_i^{(\mu)}\mid s_{-i}^{(\mu)}\bigr),$$

with estimator

$$\hat\theta_L^{PL}=\arg\max_{\theta\in\Theta} PLL_L(\theta),$$

where $s_{-i}$ denotes all components except $s_i$ (Mozeika et al., 2014). In undirected graphical models this becomes a product over nodewise conditionals, and by the local Markov property each conditional can be reduced to dependence on the Markov blanket rather than the full complement of variables (Pensar et al., 2014).

The main computational advantage is that the conditional distribution

$$P_\theta(s_i\mid s_{-i}) = \frac{P_\theta(s_i,s_{-i})}{\sum_{s_i'}P_\theta(s_i',s_{-i})}$$

does not involve the intractable normalizing constant of the joint model (Mozeika et al., 2014). This feature explains the central role of pseudo-likelihood in Gibbs random fields, Ising and Potts models, and Markov random fields, where direct likelihood maximization typically requires summation over exponentially many states or repeated MCMC evaluation of partition functions (Dikmen, 2015, Bouranis et al., 2017).

An information-theoretic interpretation sharpens this computational picture. Using the empirical distribution $P_L$, Mozeika, Dikmen, and Piili derive a relation between pseudo-likelihood and likelihood through Kullback–Leibler divergences, showing that maximizing pseudo-likelihood lower-bounds the true log-likelihood up to an additive constant (Mozeika et al., 2014). This does not make pseudo-likelihood identical to likelihood, but it explains why pseudo-likelihood can retain substantial inferential information even after replacing global normalization by empirical marginals.

2. Consistency, asymptotics, and regularity conditions

A general weak-consistency theorem for maximum pseudo-likelihood requires four conditions: identifiability of the conditionals, compactness of the parameter space, continuity of the pseudo-log-likelihood contribution $q(s,\theta)=\sum_i \log P_\theta(s_i\mid s_{-i})$, and an integrable envelope dominating $q$ uniformly in $\theta$ (Mozeika et al., 2014). Under these assumptions, the sample pseudo-likelihood $$PLL_L(\theta) = \frac{1}{L}\sum_{\mu=1}^L \sum_{i=1}^N \log P_\theta\!\bigl(s_i^{(\mu)}\mid s_{-i}^{(\mu)}\bigr),$$0 converges uniformly in probability to its population counterpart $$PLL_L(\theta) = \frac{1}{L}\sum_{\mu=1}^L \sum_{i=1}^N \log P_\theta\!\bigl(s_i^{(\mu)}\mid s_{-i}^{(\mu)}\bigr),$$1, and the maximizer $$PLL_L(\theta) = \frac{1}{L}\sum_{\mu=1}^L \sum_{i=1}^N \log P_\theta\!\bigl(s_i^{(\mu)}\mid s_{-i}^{(\mu)}\bigr),$$2 converges in probability to the true parameter $$PLL_L(\theta) = \frac{1}{L}\sum_{\mu=1}^L \sum_{i=1}^N \log P_\theta\!\bigl(s_i^{(\mu)}\mid s_{-i}^{(\mu)}\bigr),$$3 (Mozeika et al., 2014). Strong consistency can be obtained by replacing the uniform law of large numbers in probability with almost-sure uniform convergence, for example through Glivenko–Cantelli conditions (Mozeika et al., 2014).

The same paper notes that although the formal development is given for discrete variables, the arguments carry over verbatim to continuous models under the usual measurability and integrability conditions (Mozeika et al., 2014). This broader scope is reflected in later work on copula semiparametric models, where pseudo-likelihood is built from residual-based pseudo-observations rather than directly observed latent errors (Omelka et al., 2019).

In that copula setting, Omelka, Hudecová, and Neumeyer prove that, under marginal smoothness assumptions, tail control for the residual error densities, copula-family identifiability, twice continuous differentiability, and non-singular Fisher information, the estimator based on estimated residuals is asymptotically equivalent to the oracle estimator based on unobserved errors: $$PLL_L(\theta) = \frac{1}{L}\sum_{\mu=1}^L \sum_{i=1}^N \log P_\theta\!\bigl(s_i^{(\mu)}\mid s_{-i}^{(\mu)}\bigr),$$4 and therefore inherits the oracle asymptotic normality (Omelka et al., 2019). The same analysis is also a reminder of pseudo-likelihood’s sensitivity to regularity failures: when residual densities are discontinuous at a boundary and the copula score is unbounded, the maximum pseudo-likelihood estimator can exhibit large bias and variance, whereas moment-based estimators may remain stable (Omelka et al., 2019).

3. Pseudo-likelihood in graphical and inverse-covariance models

In Gaussian graphical models, pseudo-likelihood has become a central alternative to full Gaussian likelihood for sparse precision-matrix estimation. Roycraft and Rajaratnam study two convex pseudo-likelihood objectives on the cone of symmetric matrices with prescribed zeros and positive diagonal: concord,

$$PLL_L(\theta) = \frac{1}{L}\sum_{\mu=1}^L \sum_{i=1}^N \log P_\theta\!\bigl(s_i^{(\mu)}\mid s_{-i}^{(\mu)}\bigr),$$5

and conspace,

$$PLL_L(\theta) = \frac{1}{L}\sum_{\mu=1}^L \sum_{i=1}^N \log P_\theta\!\bigl(s_i^{(\mu)}\mid s_{-i}^{(\mu)}\bigr),$$6

both convex though not necessarily strictly convex on the open cone of admissible matrices (Roycraft et al., 2023).

Their main structural result is an existence-and-uniqueness criterion. For $$PLL_L(\theta) = \frac{1}{L}\sum_{\mu=1}^L \sum_{i=1}^N \log P_\theta\!\bigl(s_i^{(\mu)}\mid s_{-i}^{(\mu)}\bigr),$$7, a unique minimizer of either objective exists if and only if

$$PLL_L(\theta) = \frac{1}{L}\sum_{\mu=1}^L \sum_{i=1}^N \log P_\theta\!\bigl(s_i^{(\mu)}\mid s_{-i}^{(\mu)}\bigr),$$8

as linear subspaces; uniqueness is guaranteed in particular when both spaces reduce to $$PLL_L(\theta) = \frac{1}{L}\sum_{\mu=1}^L \sum_{i=1}^N \log P_\theta\!\bigl(s_i^{(\mu)}\mid s_{-i}^{(\mu)}\bigr),$$9 (Roycraft et al., 2023). This gives a precise answer to the long-open minimum-sample-size question for these pseudo-likelihood estimators. Defining the weak pseudo-likelihood rank $$\hat\theta_L^{PL}=\arg\max_{\theta\in\Theta} PLL_L(\theta),$$0 and the weak Gaussian rank $$\hat\theta_L^{PL}=\arg\max_{\theta\in\Theta} PLL_L(\theta),$$1, they obtain the sandwich

$$\hat\theta_L^{PL}=\arg\max_{\theta\in\Theta} PLL_L(\theta),$$2

so the pseudo-likelihood threshold never improves on the Gaussian MLE threshold on the minimal-$$\hat\theta_L^{PL}=\arg\max_{\theta\in\Theta} PLL_L(\theta),$$3 axis, although for trees, chordal graphs, grids, complete graphs, and other commonly studied classes the thresholds coincide exactly (Roycraft et al., 2023).

A second line of development modifies pseudo-likelihood to enforce strict convexity. PseudoNet estimates the inverse covariance matrix by

$$\hat\theta_L^{PL}=\arg\max_{\theta\in\Theta} PLL_L(\theta),$$4

and can be interpreted as a generalization of the Gaussian likelihood, CONCORD, and SPACE under appropriate choices of linear operators $$\hat\theta_L^{PL}=\arg\max_{\theta\in\Theta} PLL_L(\theta),$$5 (Ali et al., 2016). The additional ridge term yields uniqueness, strong convexity, linear convergence of the optimization algorithm, and prevents saturation, in contrast to pure lasso-based pseudolikelihood methods that can pick at most $$\hat\theta_L^{PL}=\arg\max_{\theta\in\Theta} PLL_L(\theta),$$6 nonzeros while PseudoNet can, for generic data matrices, pick up to $$\hat\theta_L^{PL}=\arg\max_{\theta\in\Theta} PLL_L(\theta),$$7 nonzeros (Ali et al., 2016).

For sparse binary pairwise Markov networks, the $$\hat\theta_L^{PL}=\arg\max_{\theta\in\Theta} PLL_L(\theta),$$8-regularized pseudo-likelihood can be reformulated exactly as a single sparse multiple logistic-regression problem, which allows coordinate descent and strong screening rules from generalized linear models to be applied without approximation loss (Geng et al., 2017). This reformulation yields substantial speedup and improved stability over node-wise logistic regression on unbalanced high-dimensional data when the regularization parameter is small (Geng et al., 2017).

4. Bayesian structure learning and evidence approximation

Pseudo-likelihood also supports Bayesian scoring rules for graph structure learning when full marginal likelihoods are unavailable. For discrete Markov networks, Pensar, Nyman, Corander, and Koski introduce a Bayesian marginal pseudo-likelihood obtained by integrating local conditional multinomial parameters against independent Dirichlet priors. Under global and local parameter independence, the score factorizes into closed-form Gamma-ratio terms, each depending only on local contingency counts and prior hyperparameters (Pensar et al., 2014). Asymptotically, the log-score is equivalent to the pseudo-information criterion, which yields consistency of the estimated graph under faithfulness to a true Markov network (Pensar et al., 2014). A notable feature is automatic regularization through parameter integration rather than external tuning of an $$\hat\theta_L^{PL}=\arg\max_{\theta\in\Theta} PLL_L(\theta),$$9 penalty (Pensar et al., 2014).

For Gaussian graphical models, Leppä-aho, Pensar, and Corander develop a fractional marginal pseudo-likelihood that combines Besag’s pseudo-likelihood with a fractional Bayes factor construction. Using an improper reference prior on the precision matrix and a fractional update with $s_{-i}$0, they derive closed-form local factors for each Markov blanket and obtain a tuning-free scoring criterion for arbitrary, including non-decomposable, graphs (Leppä-aho et al., 2016). Their theoretical result states that for each node the true Markov blanket uniquely maximizes the local score asymptotically, implying consistent recovery of the global undirected graph (Leppä-aho et al., 2016).

In exponential random graph models, standard pseudo-likelihood is computationally attractive but can be inferentially unreliable. Bouranis, Friel, and Maire therefore propose an adjusted pseudo-likelihood

$s_{-i}$1

constructed to match the true likelihood in mode, curvature at the mode, and magnitude at the mode (Bouranis et al., 2017). The mode adjustment aligns $s_{-i}$2 with $s_{-i}$3, the curvature adjustment uses $s_{-i}$4 from Cholesky factors of the two Hessians, and the magnitude adjustment rescales by an estimate of the intractable normalizing constant obtained through path sampling or annealed importance sampling (Bouranis et al., 2017). Empirically, the adjusted version produces Bayes factors close to gold-standard methods at much lower computational cost, whereas unadjusted pseudo-likelihood can be seriously misleading (Bouranis et al., 2017).

5. Optimization, decomposition, and distributed computation

A persistent reason for pseudo-likelihood’s longevity is algorithmic tractability. In Gaussian graphical models, concord and conspace are convex programs with linear equality constraints, and can be handled by second-order methods, coordinate descent, or ADMM (Roycraft et al., 2023). PseudoNet separates the objective into a smooth part with Lipschitz gradient and a nonsmooth $s_{-i}$5 part with soft-thresholding proximal map, so accelerated or nonaccelerated proximal-gradient iterations converge geometrically to the unique minimizer when $s_{-i}$6, with each iteration costing $s_{-i}$7 operations (Ali et al., 2016). Sequential screening rules based on derivative bounds can discard 80–90% of variables when solving across a regularization grid (Ali et al., 2016).

Pseudo-likelihood is equally useful when computation itself must be decentralized. In distributed parameter estimation for pairwise exponential-family Markov random fields, the global pseudo-likelihood decomposes into nodewise local objectives

$s_{-i}$8

and each node computes a local M-estimator using only samples restricted to its neighborhood (Liu et al., 2012). The local estimators can then be combined by linear consensus or max consensus. Under standard regularity, these combined estimators are consistent and asymptotically normal, and simple inverse-variance weighting or max-voting based on local Fisher information can be statistically competitive with joint optimization while retaining low communication cost and anytime behavior (Liu et al., 2012).

A related computational theme is the replacement of difficult derivatives by operator-based linear algebra. In efficient pseudo-likelihood estimation for dynamic discrete games, the Jacobian-free algorithm of Fukasawa avoids explicit Jacobians of equilibrium constraints by computing Jacobian-vector products through central finite differences and solving the resulting linear systems with GMRES (Fukasawa, 2024). The estimator remains the EPL estimator of Dearing and Blevins, but the implementation removes a major source of coding burden and coding error in complicated models (Fukasawa, 2024).

A further computational critique concerns regularization itself. Buchweitz, Frazier, and Martin argue that directly penalizing a pseudo-likelihood is often suboptimal because the pseudo-score may mis-specify curvature and tail behavior. Their two-stage regularization instead approximates the sampling distribution of the first-stage pseudo-MLE and regularizes that stage-two Gaussian approximation, recovering ordinary MAP estimation when the original score happens to be the true likelihood (Buchweitz et al., 2020).

6. Applications, limitations, and recent extensions

Pseudo-likelihood methods are used across a wide range of structured-data problems. In community detection for sparse stochastic block models, Amini, Chen, Bickel, and Levina replace the full network likelihood by a mixture model for block-sums computed relative to a working partition, producing unconditional and degree-conditioned pseudo-likelihood EM algorithms with computational cost $s_{-i}$9 and consistency under a mild condition on the starting labeling in the two-community case (Amini et al., 2012). The same paradigm has been extended to weighted networks with Gaussian edge weights, where a pseudo-likelihood built from block-sum Gaussian mixtures is consistent in a homogeneous weighted SBM and performs well in simulations and fMRI data (Cerqueira et al., 2023).

The method has also influenced sequence modeling and mixture learning. In Gaussian-process sequence labeling, pseudo-likelihood permits conditioning each label on the input sequence and a chosen dependency set of neighboring labels, making long-range output interactions computationally feasible under a variational Gaussian approximation (Srijith et al., 2014). For mixtures of Ising models, the mixture pseudo-likelihood replaces each joint component density by a product of single-site conditionals, leading to an EM-style coordinate-ascent procedure that can recover distinct components in synthetic data and in Potts-model applications to protein contact prediction (Dikmen, 2015).

Recent work broadens the term beyond classical graphical-model estimation. Pseudo-Likelihood Inference in simulation-based inference defines

$s_i$0

where $s_i$1 is an integral probability metric such as MMD or $s_i$2-Wasserstein distance, and adapts the bandwidth through information-theoretic trust regions (Gruner et al., 2023). This yields gradient-based training of neural posteriors without summary statistics and with support for multiple observations, and it improves on SNPE when more data are available, especially on stochastic simulations and multi-modal posterior landscapes (Gruner et al., 2023). In energy-based Ising models, maximizing pseudo-likelihood has been shown to induce zero-temperature dynamics

$s_i$3

so that a pseudo-likelihood-trained network behaves as an associative memory with large capacity and meaningful generalization even for asymmetric couplings (D'Amico et al., 7 Jul 2025).

Several misconceptions recur in the literature. Pseudo-likelihood is not, in general, the true joint likelihood; sequence-labeling work notes that estimates may therefore be biased relative to fully normalized models when global consistency matters (Srijith et al., 2014). It is also not uniformly more statistically efficient than exact likelihood; for Markov networks, Pensar et al. explicitly note that pseudo-likelihood is in general less statistically efficient than maximum likelihood and can overestimate edges (Pensar et al., 2014). Nor does it automatically improve existence properties over full Gaussian likelihood; in sparse Gaussian graphical models, Roycraft and Rajaratnam show that pseudo-likelihood never beats the Gaussian MLE on the minimal-sample-size axis, even though the thresholds coincide for many graph classes of practical interest (Roycraft et al., 2023). These caveats help explain a broader pattern: pseudo-likelihood is best understood not as a universal replacement for likelihood, but as a family of tractable surrogates whose success depends on the geometry of local conditionals, the regularity of the model class, and the fidelity with which local normalization preserves the inferential features of the original joint model.