---
title: Spatial Autoregressive Variational Distributions
url: https://www.emergentmind.com/topics/spatially-autoregressive-variational-distributions
type: topic
---

# Spatial Autoregressive Variational Distributions

Searching arXiv for the cited papers to ground the article in current literature.
Spatially autoregressive variational distributions are variational families in which the approximate posterior or latent distribution is organized so that dependence follows a spatial, graph, or locality-preserving structure rather than an exchangeable or fully factorized one. In recent work, this idea appears in several distinct but related forms: Gaussian variational posteriors whose precision is the sum of restricted-likelihood and spatial precision terms in ICAR models, checkerboard-ordered local factorizations for object catalogs in astronomical images, and variable-length autoregressive latent chains that retrieve spatially grounded information for control [2604.07635][2510.03074][2607.08724]. Related formulations use matrix-variate normal latents to encode spatial structure explicitly, variational autoencoders trained on CAR samples to emulate spatial random effects, and variational autoregressive networks for lattice systems [1705.06821][2503.14710][2409.20029].

## 1. Conceptual foundations

The term “spatially autoregressive” does not identify a single variational family. In graph-based spatial statistics, it refers to variational distributions whose covariance or precision inherits adjacency structure from a spatial weights matrix or graph Laplacian. In Gaussian ICAR/CAR models, one starts from an adjacency matrix \(W\), a degree matrix \(D=\mathrm{diag}(d_1,\dots,d_n)\), and a graph Laplacian \(Q=D-W\). Proper CAR uses a positive-definite precision, often of the form \(\tau Q_\rho\) with \(Q_\rho=D-\rho W\), whereas ICAR takes \(\rho=1\), so \(Q=D-W\) is singular for connected graphs because \(Q1=0\) and \(\mathrm{rank}(Q)=n-1\) [2604.07635].

In spatial detection and other neural settings, “autoregressive” instead refers to an ordering of spatial latent variables or tiles. In astronomical detection, the latent space is partitioned by a \(K\)-color checkerboard pattern, and the variational family factorizes over color classes,
\[
q(z\mid x)=\prod_{k=1}^K q(z_k\mid z_{<k},x),
\]
with within-rank conditional independence and local conditioning chosen to mirror the posterior Markov blanket induced by finite object radius and additive image formation [2510.03074].

In continuous-control reasoning, spatial organization can be present without a fixed grid factorization. Latent Memory Palace formulates a variable-length discrete autoregressive latent sequence \(z_{1:T(z)}\), terminated by EOS, and describes “spatial organization” as a learned space navigated sequentially while the encoder and decoder cross-attend to spatially structured observation embeddings from image encoders. The latent prior itself is sequential rather than grid-indexed, but each latent step can retrieve and focus on spatially localized cues through attention [2607.08724].

A further generalization appears in variational autoregressive networks for statistical mechanics, where a spatial configuration \(x=(x_1,\dots,x_N)\) is ordered into a sequence and modeled as
\[
q_\theta(x)=\prod_{i=1}^N q_\theta(x_i\mid x_{<i}).
\]
Here the “spatial” content is controlled by the ordering, with raster, zigzag, spiral, or space-filling curves such as Hilbert or Z-order curves used to preserve locality on lattices [2409.20029].

## 2. Graph-precision formulations in spatial statistics

The most explicit spatially autoregressive variational distributions arise when the variational covariance is tied directly to a spatial precision. In VREML for Gaussian ICAR models, the variational family is a Gaussian on the constrained subspace \(\mathcal E=\{v:1^\top v=0\}\),
\[
q(u)=\mathcal N_{\mathcal E}(m,S),
\]
and the ELBO on the restricted likelihood simplifies to
\[
\begin{aligned}
\mathcal L_V(m,S,\tau_y,\tau_u)
&=\ \tfrac{n-p}{2}\log \tau_y
-\tfrac{\tau_y}{2}\Big[(y-m)^\top K(y-m)+\operatorname{tr}(K S)\Big]\\
&\quad+\tfrac{n-r}{2}\log \tau_u
-\tfrac{\tau_u}{2}\Big[m^\top Q m + \operatorname{tr}(Q S)\Big]
+\tfrac{1}{2}\log|S|_+ + \text{const.}
\end{aligned}
\]
The defining property is that the optimal variational precision is
\[
S^{-1}=\tau_y K+\tau_u Q,
\]
so the spatial autoregressive precision \(Q\) appears both in the ELBO and in the variational covariance itself. Under Gaussian \(y\mid \beta,u\) and an ICAR prior on \(u\), the conditional posterior \(p(u\mid y,\tau_y,\tau_u)\) is Gaussian on \(\mathcal E\) with the same precision and mean \(m=S\tau_y K y\); consequently, the variational family is exact and the ELBO is tight [2604.07635].

Spatial error models provide a related but distinct construction. In the SEM,
\[
y=X\beta+v,\qquad v=\rho Wv+\varepsilon,
\]
with \(A=I_n-\rho W\), the induced covariance is \(\sigma_y^2(A^\top A)^{-1}\) and the spatial precision is \(M_y=A^\top A\). Variational Bayes with missing responses uses Gaussian variational approximations with factor covariance,
\[
q_\lambda(\theta,y_u)=\mathcal N\big((\theta^\top,y_u^\top)^\top;\ \mu,\ BB^\top + D^2\big),
\]
or a hybrid factorization
\[
q_\lambda(\theta,y_u)=p(y_u\mid O,\theta)\,q_{\lambda^0}(\theta).
\]
Although the variational covariance is low-rank-plus-diagonal rather than explicitly sparse-precision, spatial dependence enters the optimization through \(\log|M_y|\), \((y-X\beta)^\top M_y(y-X\beta)\), and the exact conditional
\[
p(y_u\mid \phi,y_o)=\mathcal N\!\Big(X_u\beta - M_{y,uu}^{-1}M_{y,uo}(y_o - X_o\beta),\ \sigma_y^2 M_{y,uu}^{-1}\Big),
\]
which is used by HVB under MAR [2406.08685].

The same SEM backbone has been extended to non-Gaussian settings by combining Student’s \(t\)-distributed errors with the Yeo–Johnson transformation. The paper introduces SEM-t, YJ-SEM-Gau, and YJ-SEM-t, and develops VB and HVB estimators based on a Gaussian variational approximation with factor covariance and reparameterized constrained parameters \(\omega'=\log(\sigma_e^2)\), \(\rho'=\log(1+\rho)-\log(1-\rho)\), \(\nu'=\log(\nu-3)\), and \(\gamma'=\log(\gamma)-\log(2-\gamma)\). In these models, the spatial autoregressive structure continues to operate through \(A=I_n-\rho W\), the determinant term \(\log|A^\top A|\), and quadratic forms such as \(r^\top A^\top A r\) or \(r^\top A^\top \Sigma_\tau^{-1} A r\) [2505.23070].

A further generalization appears in unrestricted panel SAR models. The spatial lag is represented by a learned \(N\times N\) spillover matrix \(\Lambda\) in
\[
y_t=\Lambda y_t + \widetilde X_t \widetilde\beta + u_t,
\]
and variational Bayes is carried out row-by-row in two stages. Stage 1 estimates \(Y_{/i}=X\Upsilon_i+E_i\) with Gaussian \(q(\gamma)\) for \(\gamma=\mathrm{vec}(\Upsilon_i)\) and a learned precision matrix \(\Omega\) for the first-stage errors. Stage 2 estimates
\[
y_i=\widehat Y_{/i}(\Lambda_{i\cdot})' + X_{\cdot,i}\widetilde\beta_i + u_i
\]
with Gaussian \(q(\theta)\) for \(\theta=[(\Lambda_{i\cdot})'\ \widetilde\beta_i]\). Dirichlet–Laplace priors induce sparsity in both \(\Lambda\) and \(\Omega\), yielding closed-form variational updates for Gaussian, Gamma, inverse Gaussian, and generalized inverse Gaussian factors [2205.15420].

## 3. Matrix-normal and autoencoded spatial latents

Spatial variational auto-encoding via matrix-variate normal distributions provides a different route to spatially structured variational inference. Instead of vector latents, the latent variable is a feature map \(Z\in\mathbb R^{m\times n}\) with
\[
Z\sim MN(M,U,V),
\qquad
\mathrm{vec}(Z)\sim N(\mathrm{vec}(M),V\otimes U).
\]
The Kronecker covariance \(V\otimes U\) imposes a separable spatial structure on the latent feature map, with \(U\) controlling covariance across rows and \(V\) across columns. Within the VAE ELBO,
\[
\mathrm{ELBO}(x)=\mathbb E_{q(Z\mid x)}[\log p(x\mid Z)]-\mathrm{KL}(q(Z\mid x)\|p(Z)),
\]
the KL against a standard normal prior reduces to a Gaussian KL with \(\Sigma_q=V\otimes U\), and the paper develops a low-rank mean parameterization \(M_k=\mu_k\nu_k^\top\) that reduces encoder outputs from \((d^2+2d)N\) to \(4dN\) per collection of \(N\) latent maps of size \(d\times d\). The paper explicitly notes that it does not use AR(1) Toeplitz or GMRF/SAR precision structures, but also states that such autoregressive covariance choices are natural extensions of the MVN framework [1705.06821].

Variational autoencoded multivariate spatial Fay–Herriot models move the spatial autoregressive burden into a learned generator. In the separable MCAR case,
\[
\mathrm{vec}(\boldsymbol{\phi})\sim
\mathcal N\!\big(\mathbf 0,\ (\boldsymbol D-\rho \boldsymbol W)^{-1}\otimes \boldsymbol\Sigma\big),
\]
which the paper rewrites as
\[
\boldsymbol{\psi}\sim MN_{N\times K}\big(\mathbf 0,\ (\boldsymbol D-\rho \boldsymbol W)^{-1},\ \boldsymbol I\big),
\qquad
\boldsymbol{\phi}=\boldsymbol{\psi}\boldsymbol L^\top.
\]
A \(\beta\)-VAE is then trained on samples from the CAR prior, with latent standard normal \(z\) and decoder \(\mathrm{decoder}_\zeta(z)\). After training, the decoder is fixed and reused as a prior emulator for spatial random effects on the same geography. For VSMS-FH, the implementation uses \(\tilde{\boldsymbol\phi}=\sigma\,\tilde{\boldsymbol\psi}\); for VGMS-FH, two decoder draws are combined through the GMCAR structure
\[
\tilde{\boldsymbol\phi}_2=\sigma_2\,\mathrm{decoder}_\zeta(\boldsymbol z_2),
\qquad
\tilde{\boldsymbol\phi}_1=\boldsymbol A\,\tilde{\boldsymbol\phi}_2+\sigma_1\,\mathrm{decoder}_\zeta(\boldsymbol z_1),
\]
with \(\boldsymbol A=\eta_0\boldsymbol I+\eta_1\boldsymbol W\). Spatial dependence is therefore not imposed by an explicit variational precision during downstream inference; it is learned from samples generated with \(\boldsymbol Q(\rho)=\boldsymbol D-\rho \boldsymbol W\), using \(\rho\sim\mathrm{Unif}(0,1)\) during VAE training [2503.14710].

These matrix-normal and autoencoded constructions shift attention from sparse-precision exactness to amortized spatial structure. One uses explicit covariance algebra in latent feature maps, and the other uses a learned decoder as a reusable generator of spatial effects.

## 4. Checkerboard and memory-palace autoregression

In high-dimensional perception and control, spatially autoregressive variational distributions are often built around locality and conditional independence rather than sparse Gaussian precision. For object detection in astronomical images, the latent catalog is tiled spatially and split into tiles of interest \(z_\ell\). Because object contributions are local—\(\xi_{h,w}(u_s,v_s)=0\) when \(\|[h,w]^\top-u_s\|_\infty>R\)—the posterior has a local Markov structure. The variational family uses a \(K\)-color checkerboard partition with rank
\[
\Psi(h,w)=((h-1)\bmod \sqrt K)\cdot \sqrt K + ((w-1)\bmod \sqrt K),
\]
and factorizes as
\[
q(z\mid x)=\prod_{k=1}^K \prod_{\ell\in C_k}
q(z_\ell\mid z_{<k}\cap z_\ell^{r_{\mathcal N}},\,x_{(\ell-1)T+1}^{r_{\mathcal X}}).
\]
For a posterior locality radius \(R\) and tile size \(T\), the minimal number of colors is \(\sqrt K=\lceil 2R/T+1\rceil\). Within each tile, an auxiliary permutation \(\sigma\) gives
\[
q(z_k\mid z_{<k},x)
=
\mathbb E_{\sigma\sim q(\sigma)}
\Big[\prod_{i=1}^M q(z_k^{[i]}\mid z_k^{[<i]},z_{<k},x)\Big].
\]
This yields a variational graph that mirrors the posterior conditional independencies while preserving GPU parallelism within each color class [2510.03074].

Latent Memory Palace adopts a variable-length discrete autoregressive latent sequence rather than a tiled field. With observation \(o\), action chunk \(a\), latent vocabulary \(V\cup\{\mathrm{EOS}\}\), and termination time \(T(z)=\min\{t\ge 1:z_t=\mathrm{EOS}\}\), the model uses
\[
p_\theta(z\mid o)=\prod_{t=1}^{T(z)} p_\theta(z_t\mid z_{<t},o),
\qquad
q_\theta(z\mid o,a)=\prod_{t=1}^{T(z)} q_\theta(z_t\mid z_{<t},o,a),
\]
with decoder \(p_\phi(a\mid o,z)=p_\phi(a\mid o,z_{1:T(z)})\). The latent sequence is interpreted as an iterative retrieval path in a learned “memory palace.” The paper emphasizes that this is not a PixelCNN-style factorization over a fixed image grid and not an explicit geometric graph. Instead, latent steps form a route through a learned space, and spatial grounding is supplied by cross-attention to visual observation tokens and positional embeddings. The resulting latent path is therefore sequential in its prior factorization but spatial in what it retrieves [2607.08724].

These two examples illustrate a central distinction. In the checkerboard family, spatial autoregression is made explicit by coloring and neighborhood conditioning. In the memory-palace family, spatial structure is implicit in the retrieval mechanism and the observation embeddings, while autoregression governs the order and length of latent reasoning.

## 5. Optimization, training criteria, and inference mechanics

The optimization of spatially autoregressive variational distributions is as diverse as their parameterizations. In VREML, the ELBO admits closed-form coordinate-ascent updates:
\[
S^{-1}=\tau_y K+\tau_u Q,\qquad
m=S\tau_y K y,
\]
\[
\tau_y=\frac{n-p}{(y-m)^\top K (y-m)+\operatorname{tr}(K S)},
\qquad
\tau_u=\frac{n-r}{m^\top Q m+\operatorname{tr}(Q S)}.
\]
The paper proves that each block is strictly concave, the ELBO is monotone nondecreasing under block updates, and any accumulation point is a stationary point [2604.07635].

Astronomical SAVD is trained not by an ELBO on observed data but by Neural Posterior Estimation under the expected forward KL. The objective is
\[
\mathcal L(\phi)=-\mathbb E_{p(z)p(x\mid z)}[\log q_\phi(z\mid x)],
\]
with unbiased gradient
\[
\nabla \mathcal L(\phi)=-\mathbb E_{p(z)p(x\mid z)}[\nabla \log q_\phi(z\mid x)].
\]
The paper stresses the contrast between forward and reverse KL: forward KL is mode-covering and better aligned with calibration for ambiguous blends, while the autoregressive tiling ensures that the variational graph matches the posterior’s locality [2510.03074].

LMP uses neither continuous reparameterization nor Gumbel-softmax. Because discrete autoregressive sampling blocks direct reparameterization, the ELBO is optimized by latent-space reinforcement learning. The KL term decomposes stepwise,
\[
D_{\mathrm{KL}}(q_\theta(z\mid o,a)\|p_\theta(z\mid o))
=
\sum_{t=1}^\infty
\widetilde{\mathbb E}_{z_{\le t}\sim q_\theta(\cdot\mid o,a)}
\big[
\log q_\theta(z_t\mid z_{<t},o,a)-\log p_\theta(z_t\mid z_{<t},o)
\big],
\]
which yields per-step rewards. Stable updates use a PPO-style clipped surrogate with length-aware clipping, prior mixing, and per-step free-nats. Adaptive halting is induced by a decoder variance schedule \(\sigma(T)=\gamma^T\sigma_0\), so extra latent steps are useful only when they improve the reconstruction mean enough to offset tighter variance [2607.08724].

For variational autoregressive networks in statistical mechanics, the objective is variational free energy rather than a standard data ELBO. With target Boltzmann distribution \(p(x)\propto \exp(-\beta E(x))\), minimizing \(KL(q_\theta\|p)\) is equivalent to minimizing
\[
F(q_\theta)=\mathbb E_{q_\theta}[E(x)]-\frac{1}{\beta}H(q_\theta).
\]
The natural-gradient update uses the Fisher information
\[
F(\theta)=\mathbb E_{q_\theta}[s(x)s(x)^\top],\qquad s(x)=\nabla_\theta \log q_\theta(x),
\]
and, in low-rank sample space, becomes
\[
\Delta \theta
=
-\alpha O^\top(OO^\top+\xi I)^{-1}R.
\]
This changes the dominant inversion from parameter dimension \(P\) to batch size \(B\), giving complexity cubic in batch size rather than in model parameters [2409.20029].

Missing-data spatial models often require hybridization. In SEMs under MAR or MNAR, HVB factorizes as \(q(\theta,y_u)=p(y_u\mid O,\theta)q(\theta)\) or its MNAR analogue and draws \(y_u\) by exact Gaussian conditionals, block Gibbs, or Metropolis–Hastings proposals embedded inside the variational loop. This approach is used both in Gaussian SEMs and in non-Gaussian YJ-SEM and \(t\)-SEM variants [2406.08685][2505.23070]. In unrestricted panel SAR, the corresponding hybridization is algebraic rather than stochastic: two-stage VB decomposes the spatial simultaneous system into first-stage instrument construction and second-stage row-wise estimation, with Dirichlet–Laplace shrinkage making the updates closed form [2205.15420].

## 6. Empirical behavior, misconceptions, and open questions

Empirical results show that the benefit of spatially autoregressive variational structure depends strongly on whether the underlying task is governed by local dependence, adjacency-constrained smoothness, or multimodal spatial retrieval. In Xenium breast cancer spatial transcriptomics, VREML reports RMSE \(0.5752\) and MAE \(0.4045\), compared with Exact REML \(0.5984\), \(0.4196\), MLE \(0.5983\), \(0.4196\), and INLA \(0.5982\), \(0.4195\). The same paper states that VREML is “exact under Gaussian ICAR settings,” so the empirical advantage is computational and numerical rather than an approximation gap in \(u\) [2604.07635].

In astronomical detection, the checkerboard autoregressive structure is tied directly to calibration. On held-out synthetic SDSS-like data, log-likelihood of ground-truth catalogs improves from \(-227{,}337\) for \(K=1\) to \(-210{,}833\) for \(K=4\), and F1 improves from approximately \(0.910\) to \(0.923\). In the crowded M2 setting, log-likelihood improves from \(-4{,}216{,}070\) to \(-3{,}454{,}335\), and F1 improves from approximately \(0.571\) to \(0.589\). The paper repeatedly attributes these gains to reduced border artifacts and better calibration near tile boundaries [2510.03074].

In control, LMP-\(\pi\) and LMP-\(\texttt{tok}\) connect spatial retrieval to adaptive test-time compute. On DROID zero-shot tasks, LMP-\(\pi\) reports block-bowl \(0.65\) versus \(0.40\) for DP and marker-mug \(0.55\) versus \(0.25\). On finetuned DROID tasks, it reports peg-hole \(0.70\) versus \(0.30\), and clean-table \(1/3\), \(2/3\), \(3/3\) values of \(0.95\) versus \(0.80\), \(0.95\) versus \(0.55\), and \(0.55\) versus \(0.25\). On LIBERO-90 multitask, average performance is \(0.933\pm 0.003\) versus \(0.909\pm 0.004\), and bottom-10 tasks are \(0.645\pm 0.023\) versus \(0.463\pm 0.011\). The paper further reports a negative correlation \(r=-0.518\) between average latent steps and KNN action variance, and interprets this as linking adaptive compute to irreducible uncertainty [2607.08724].

Autoencoded spatial priors in Fay–Herriot models are motivated primarily by scalability. In Missouri simulations, SMS-FH requires about \(12.72\) hours per simulation and VSMS-FH about \(7.88\) minutes, while GMS-FH requires about \(3.14\) hours and VGMS-FH about \(6.47\) minutes. In California, the paper states that GMS-FH would require more than \(45\) days for one simulation, whereas VGMS-FH requires about \(2.80\) hours and VSMS-FH about \(3.51\) hours per simulation on GPU. The empirical pattern reported is that the VAE-based models preserve most of the spatial benefit of MCAR or GMCAR structure while making tract-scale multivariate estimation computationally feasible [2503.14710].

For SEMs with missing data, the main practical lesson is that variational structure alone is not always sufficient. Under MAR with \(n=10{,}000\) and \(75\%\) missing, JVB gives \(\hat{\sigma}_y^2\approx 2.15\) and \(\hat{\rho}\approx 0.08\), whereas HVB-G gives \(\hat{\sigma}_y^2\approx 1.00\) and \(\hat{\rho}\approx 0.80\). Under MNAR with the same dimensions and missingness, HVB-AllB and HVB-3B give \(\hat{\rho}\approx 0.81\) and \(0.82\), versus JVB \(\hat{\rho}\approx 0.14\). In the non-Gaussian SEM literature, YJ-SEM-t is reported to provide the best fit for skewed heavy-tailed data, while YJ-SEM-Gau is often nearly as competitive at lower computational expense [2406.08685][2505.23070].

Two recurrent misconceptions are explicitly rejected by this literature. First, spatially autoregressive variational inference is not synonymous with gridwise factorization: LMP states that its latent autoregression is “over sequence positions that represent a navigable route in memory, not fixed spatial coordinates,” while SAVD uses color classes and local neighborhoods rather than a pixelwise raster. Second, variational inference in spatial models is not always inexact: in Gaussian ICAR restricted-likelihood estimation, the optimal Gaussian variational family coincides with the true restricted posterior [2607.08724][2510.03074][2604.07635].

Open directions are stated in several papers. LMP identifies sensitivity to RL hyperparameters, latent collapse, and the cost of larger vocabularies and longer horizons, and suggests continuous latent chains, persistent latents across an episode, and hierarchical spatial-temporal factorizations. SAVD notes exposure bias in NPE, the trade-off between larger \(K\) and sequential cost, and possible extensions through multi-scale tilings or learned neighborhoods. The SEM and panel-SAR papers emphasize scalability of blocked MCMC, weak-instrument issues, and the possibility of variational families with sparse-plus-low-rank precisions that mirror the spatial operator more directly. A plausible implication is that future work will continue to move between two poles: exact or near-exact variational families tied to spatial precision operators, and amortized neural families that learn spatial dependence from data while sacrificing direct probabilistic structure [2607.08724][2510.03074][2406.08685][2205.15420].

Source: https://www.emergentmind.com/topics/spatially-autoregressive-variational-distributions