Virtual Residual Likelihood

• Virtual Residual Likelihood reframes classical REML in Gaussian ICAR models as a variational optimization problem, enabling exact and scalable inference.
• It derives a tractable evidence lower bound (ELBO) using a Gaussian variational family, facilitating closed-form coordinate ascent updates for both mean and covariance.
• Empirical results demonstrate that VREML achieves comparable or superior accuracy and speed compared to INLA and traditional REML, especially for large spatial datasets.

Virtual (Variational) Residual Likelihood is a computational methodology for scalable inference in Gaussian spatial models, particularly those structured by intrinsic conditional autoregressive (ICAR) priors. The approach reframes the classic restricted maximum likelihood (REML) estimation as a variational optimization problem, deriving a tractable evidence lower bound (ELBO) that can be efficiently maximized via coordinate ascent. Unlike conventional REML or integrated nested Laplace approximation (INLA), virtual residual likelihood (VREML) delivers computational speedups, exactness in the Gaussian ICAR context, and closed-form updates, all while preserving the statistical properties of the restricted likelihood (Thakur, 8 Apr 2026).

1. Classical REML for Gaussian ICAR Models

Gaussian mixed models with ICAR random effects serve as the foundational setting. Given response vector $Y \in \mathbb{R}^n$, fixed-effect design $X \in \mathbb{R}^{n \times p}$, latent spatial effects $u \in \mathbb{R}^n$ (subject to the sum-to-zero constraint $1^\top u = 0$), observational precision $\tau_y > 0$, and ICAR precision $\tau_u > 0$, the model is specified as:

• $$Y|\beta, u, \tau_y \sim N(X\beta + u, \tau_y^{-1}I_n)$$,
• $$u|\tau_u \propto \tau_u^{(n - r)/2}\exp(-\frac{1}{2} \tau_u u^\top R u)$$ with $1^\top u = 0$, where $R$ is the ICAR precision matrix and $X \in \mathbb{R}^{n \times p}$0 its rank deficiency.

The restricted likelihood, integrating out $X \in \mathbb{R}^{n \times p}$1 under a flat prior, is

$X \in \mathbb{R}^{n \times p}$2

and the restricted log-likelihood is given, up to constants, as

$X \in \mathbb{R}^{n \times p}$3

where $X \in \mathbb{R}^{n \times p}$4 (Thakur, 8 Apr 2026).

2. Variational Lower Bound (ELBO) Construction

The core innovation of VREML is substituting direct maximization of $X \in \mathbb{R}^{n \times p}$5 with optimization of a variational lower bound. Introducing a variational density $X \in \mathbb{R}^{n \times p}$6, for any choice of $X \in \mathbb{R}^{n \times p}$7,

$X \in \mathbb{R}^{n \times p}$8

and by Jensen's inequality,

$X \in \mathbb{R}^{n \times p}$9

Maximizing $u \in \mathbb{R}^n$0 within an appropriate family yields the tightest lower bound; the difference $u \in \mathbb{R}^n$1 is precisely the KL divergence $u \in \mathbb{R}^n$2 (Thakur, 8 Apr 2026).

3. Gaussian Variational Family and Closed-form Objective

Choosing $u \in \mathbb{R}^n$3 to be Gaussian on the subspace $u \in \mathbb{R}^n$4, i.e., $u \in \mathbb{R}^n$5 with $u \in \mathbb{R}^n$6, $u \in \mathbb{R}^n$7, renders all required expectations in the ELBO analytically tractable:

• $u \in \mathbb{R}^n$8,
• $u \in \mathbb{R}^n$9,
• $1^\top u = 0$0, where $1^\top u = 0$1 is the pseudo-determinant on $1^\top u = 0$2.

The resulting variational lower bound, as a function of $1^\top u = 0$3, becomes

$1^\top u = 0$4

\sup_{q \in Q}\mathcal{L}V(q, \theta) = \ell{RE}(\theta),

$1^\top u = 0$5

implying that VREML attains the genuine restricted likelihood value, nullifying posterior approximation error in this setting (Thakur, 8 Apr 2026).

7. Computational Complexity and Empirical Performance

Classical REML and INLA approaches perform repeated factorization of large $1^\top u = 0$6 sparse matrices, costing $1^\top u = 0$7 per factorization for 2D lattices, with multiple such operations per parameter update. VREML, by contrast, requires a single update of the same system per iteration—reusing the factorization for both $1^\top u = 0$8 and $1^\top u = 0$9—yielding each iteration at $\tau_y > 0$0 complexity, with practical convergence reached in approximately 10–20 iterations.

Empirical benchmarks indicate that on grid sizes $\tau_y > 0$1 to $\tau_y > 0$2, VREML converges in a few seconds, while INLA or REML require substantially longer (tens of seconds to minutes). In simulation and real-data examples—including breast-cancer gene-expression data—VREML matches or slightly outperforms INLA and exact REML in mean-squared and mean absolute prediction error (MSPE/MAE) for both $\tau_y > 0$3 and $\tau_y > 0$4, with comparable variance component estimation (RMSE), but substantial gains in computational scalability (Thakur, 8 Apr 2026).

In summary, VREML substitutes the computationally intensive Gaussian restricted likelihood with a variational lower bound that is exact for Gaussian ICAR priors, yielding closed-form, monotone-convergent updates, and greatly improving scalability for large areal spatial data.