---
title: Virtual Residual Likelihood
url: https://www.emergentmind.com/topics/virtual-residual-likelihood
type: topic
---

# Virtual Residual Likelihood

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 [2604.07635].

## 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 $r$ its rank deficiency.

The restricted likelihood, integrating out $\beta$ under a flat prior, is
$$
L(\tau_y, \tau_u; Y) = \int_{\mathbb{R}^n} g(Y, u | \tau_y, \tau_u)\ du,
$$
and the restricted log-likelihood is given, up to constants, as
$$
\ell_{RE}(\theta) = -\frac{1}{2}[(n - p)\log \tau_y + \log|X^\top X| + \tau_y Y^\top P_X^\perp Y]
+ \frac{n - r}{2}\log \tau_u - \frac{1}{2}\tau_u\ \cdots,
$$
where $P_X^\perp = I - X(X^\top X)^{-1}X^\top$ [2604.07635].

## 2. Variational Lower Bound (ELBO) Construction

The core innovation of VREML is substituting direct maximization of $\ell_{RE}$ with optimization of a variational lower bound. Introducing a variational density $q(u)$, for any choice of $q$,
$$
\ell_{RE}(\tau_y, \tau_u) = \log \int g(Y, u | \tau_y, \tau_u)\ du = \log \int q(u)\frac{g(Y, u | \tau_y, \tau_u)}{q(u)} du,
$$
and by Jensen's inequality,
$$
\ell_{RE}(\tau_y, \tau_u) \geq \mathbb{E}_q[\log g(Y, u | \tau_y, \tau_u) - \log q(u)] =: \mathcal{L}_V(q; \tau_y, \tau_u).
$$
Maximizing $\mathcal{L}_V$ within an appropriate family yields the tightest lower bound; the difference $\ell_{RE} - \mathcal{L}_V$ is precisely the KL divergence $KL(q(u)\|p(u|Y, \theta))$ [2604.07635].

## 3. Gaussian Variational Family and Closed-form Objective

Choosing $q(u)$ to be Gaussian on the subspace $\mathcal{E} = \left\{ u: 1^\top u = 0 \right\}$, i.e., $q(u) = N_\mathcal{E}(\mu, \Sigma)$ with $\Sigma 1 = 0$, $1^\top \mu = 0$, renders all required expectations in the ELBO analytically tractable:
- $\mathbb{E}_q[(Y-u)^\top P_X^\perp (Y-u)] = (Y - \mu)^\top P_X^\perp (Y - \mu) + \operatorname{tr}(P_X^\perp \Sigma)$,
- $\mathbb{E}_q[u^\top R u] = \mu^\top R \mu + \operatorname{tr}(R\Sigma)$,
- $-\mathbb{E}_q[\log q(u)] = \frac{1}{2}\log|\Sigma|_+ + \frac{n - r}{2}(1 + \log 2\pi)$,
where $|\Sigma|_+$ is the pseudo-determinant on $\mathcal{E}$.

The resulting variational lower bound, as a function of $\mu, \Sigma, \tau_y, \tau_u$, becomes
$$
\mathcal{L}_V(\mu, \Sigma, \tau_y, \tau_u) = \frac{n-p}{2}\log \tau_y - \frac{\tau_y}{2}\left( (Y - \mu)^\top P_X^\perp (Y - \mu) + \operatorname{tr}(P_X^\perp\Sigma) \right)
+ \frac{n - r}{2}\log \tau_u - \frac{\tau_u}{2} \left( \mu^\top R\mu + \operatorname{tr}(R\Sigma) \right)
+ \frac{1}{2}\log|\Sigma|_+ + \mathrm{const}$ [2604.07635].

## 4. Coordinate-Ascent Updates and Algorithm

Block-wise maximization of $\mathcal{L}_V$ leads to explicit update equations:
- $\mu \leftarrow (\tau_y P_X^\perp + \tau_u R)_*^{-1} (\tau_y P_X^\perp Y)$,
- $\Sigma \leftarrow (\tau_y P_X^\perp + \tau_u R)_*^{-1}$,
- $\tau_y \leftarrow (n - p)/ \left[(Y - \mu)^\top P_X^\perp(Y - \mu) + \operatorname{tr}(P_X^\perp \Sigma)\right]$,
- $\tau_u \leftarrow (n - r)/\left[ \mu^\top R\mu + \operatorname{tr}(R\Sigma) \right]$,

where $(\cdot)_*^{-1}$ denotes the pseudo-inverse on $\mathcal{E}$. The updates proceed cyclically until convergence. The minimal form of the algorithm is as follows:

```python
Initialize τ_y⁰, τ_u⁰>0
For t=0,1,2,… until convergence:
    Σᵗ⁺¹ ← (τ_yᵗ P_X^⊥ + τ_uᵗ R)_*⁻¹
    μᵗ⁺¹ ← Σᵗ⁺¹ (τ_yᵗ P_X^⊥ Y)
    τ_yᵗ⁺¹← (n−p) / [ (Y−μᵗ⁺¹)ᵀP_X^⊥(Y−μᵗ⁺¹) + tr(P_X^⊥ Σᵗ⁺¹) ]
    τ_uᵗ⁺¹← (n−r) / [ μᵗ⁺¹ᵀR μᵗ⁺¹ + tr(R Σᵗ⁺¹) ]
```
[2604.07635]

## 5. Monotonicity and Convergence of the Variational Lower Bound

Under suitable regularity—full-rank $X$, a connected neighborhood structure $R$, and strict positivity of $\tau_y P_X^\perp + \tau_u R$ on $\mathcal{E}$—the variational lower bound $\mathcal{L}_V$ is strictly concave in each block parameter when holding others fixed, and each coordinate update increases (or preserves) $\mathcal{L}_V$. The sequence $\mathcal{L}_V^{(t)}$ generated by the coordinate-ascent updates is non-decreasing and bounded above by $\ell_{RE}$, guaranteeing convergence to a stationary point of $\mathcal{L}_V$ [2604.07635].

## 6. Exactness in the Gaussian ICAR Setting

In the case of Gaussian ICAR models, the conditional posterior $p(u|Y, \tau_y, \tau_u)$ is Gaussian on $\mathcal{E}$, coinciding exactly with the chosen variational family. Thus, the KL divergence vanishes, meaning $q(u)$ can match $p(u|Y, \tau_y, \tau_u)$ exactly:
$$
\sup_{q \in Q}\mathcal{L}_V(q, \theta) = \ell_{RE}(\theta),
$$
implying that VREML attains the genuine restricted likelihood value, nullifying posterior approximation error in this setting [2604.07635].

## 7. Computational Complexity and Empirical Performance

Classical REML and INLA approaches perform repeated factorization of large $n \times n$ sparse matrices, costing $\mathcal{O}(n^{1.5})$ 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 $\mu$ and $\Sigma$—yielding each iteration at $\mathcal{O}(n^{1.5})$ complexity, with practical convergence reached in approximately 10–20 iterations.

Empirical benchmarks indicate that on grid sizes $n=50^2$ to $200^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 $Y$ and $u$, with comparable variance component estimation (RMSE), but substantial gains in computational scalability [2604.07635].

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.

Source: https://www.emergentmind.com/topics/virtual-residual-likelihood