Spatial Autoregressive Variational Distributions

• Spatially autoregressive variational distributions are models that incorporate spatial dependence via structured covariance and autoregressive latent orderings.
• They combine graph-based precision formulations with localized factorization methods to accurately capture spatial dependencies in diverse applications.
• This approach enables scalable inference by balancing computational efficiency with the preservation of spatial correlations in latent representations.

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 (Thakur, 8 Apr 2026, Regier, 3 Oct 2025, Zhu et al., 9 Jul 2026). 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 (Wang et al., 2017, Wang et al., 18 Mar 2025, Liu et al., 2024).

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$ (Thakur, 8 Apr 2026).

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,

$D=\mathrm{diag}(d_1,\dots,d_n)$0

with within-rank conditional independence and local conditioning chosen to mirror the posterior Markov blanket induced by finite object radius and additive image formation (Regier, 3 Oct 2025).

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 $D=\mathrm{diag}(d_1,\dots,d_n)$1, 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 (Zhu et al., 9 Jul 2026).

A further generalization appears in variational autoregressive networks for statistical mechanics, where a spatial configuration $D=\mathrm{diag}(d_1,\dots,d_n)$2 is ordered into a sequence and modeled as

$D=\mathrm{diag}(d_1,\dots,d_n)$3

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 (Liu et al., 2024).

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 $D=\mathrm{diag}(d_1,\dots,d_n)$4,

$D=\mathrm{diag}(d_1,\dots,d_n)$5

and the ELBO on the restricted likelihood simplifies to

$D=\mathrm{diag}(d_1,\dots,d_n)$6

The defining property is that the optimal variational precision is

$D=\mathrm{diag}(d_1,\dots,d_n)$7

so the spatial autoregressive precision $D=\mathrm{diag}(d_1,\dots,d_n)$8 appears both in the ELBO and in the variational covariance itself. Under Gaussian $D=\mathrm{diag}(d_1,\dots,d_n)$9 and an ICAR prior on $Q=D-W$0, the conditional posterior $Q=D-W$1 is Gaussian on $Q=D-W$2 with the same precision and mean $Q=D-W$3; consequently, the variational family is exact and the ELBO is tight (Thakur, 8 Apr 2026).

Spatial error models provide a related but distinct construction. In the SEM,

$Q=D-W$4

with $Q=D-W$5, the induced covariance is $Q=D-W$6 and the spatial precision is $Q=D-W$7. Variational Bayes with missing responses uses Gaussian variational approximations with factor covariance,

$Q=D-W$8

or a hybrid factorization

$Q=D-W$9

Although the variational covariance is low-rank-plus-diagonal rather than explicitly sparse-precision, spatial dependence enters the optimization through $\tau Q_\rho$0, $\tau Q_\rho$1, and the exact conditional

$\tau Q_\rho$2

which is used by HVB under MAR (Wijayawardhana et al., 2024).

The same SEM backbone has been extended to non-Gaussian settings by combining Student’s $\tau Q_\rho$3-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 $\tau Q_\rho$4, $\tau Q_\rho$5, $\tau Q_\rho$6, and $\tau Q_\rho$7. In these models, the spatial autoregressive structure continues to operate through $\tau Q_\rho$8, the determinant term $\tau Q_\rho$9, and quadratic forms such as $Q_\rho=D-\rho W$0 or $Q_\rho=D-\rho W$1 (Wijayawardhana et al., 29 May 2025).

A further generalization appears in unrestricted panel SAR models. The spatial lag is represented by a learned $Q_\rho=D-\rho W$2 spillover matrix $Q_\rho=D-\rho W$3 in

$Q_\rho=D-\rho W$4

and variational Bayes is carried out row-by-row in two stages. Stage 1 estimates $Q_\rho=D-\rho W$5 with Gaussian $Q_\rho=D-\rho W$6 for $Q_\rho=D-\rho W$7 and a learned precision matrix $Q_\rho=D-\rho W$8 for the first-stage errors. Stage 2 estimates

$Q_\rho=D-\rho W$9

with Gaussian $\rho=1$0 for $\rho=1$1. Dirichlet–Laplace priors induce sparsity in both $\rho=1$2 and $\rho=1$3, yielding closed-form variational updates for Gaussian, Gamma, inverse Gaussian, and generalized inverse Gaussian factors (Gefang et al., 2022).

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 $\rho=1$4 with

$\rho=1$5

The Kronecker covariance $\rho=1$6 imposes a separable spatial structure on the latent feature map, with $\rho=1$7 controlling covariance across rows and $\rho=1$8 across columns. Within the VAE ELBO,

$\rho=1$9

the KL against a standard normal prior reduces to a Gaussian KL with $Q=D-W$0, and the paper develops a low-rank mean parameterization $Q=D-W$1 that reduces encoder outputs from $Q=D-W$2 to $Q=D-W$3 per collection of $Q=D-W$4 latent maps of size $Q=D-W$5. 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 (Wang et al., 2017).

Variational autoencoded multivariate spatial Fay–Herriot models move the spatial autoregressive burden into a learned generator. In the separable MCAR case,

$Q=D-W$6

which the paper rewrites as

$Q=D-W$7

A $Q=D-W$8-VAE is then trained on samples from the CAR prior, with latent standard normal $Q=D-W$9 and decoder $Q1=0$0. 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 $Q1=0$1; for VGMS-FH, two decoder draws are combined through the GMCAR structure

$Q1=0$2

with $Q1=0$3. Spatial dependence is therefore not imposed by an explicit variational precision during downstream inference; it is learned from samples generated with $Q1=0$4, using $Q1=0$5 during VAE training (Wang et al., 18 Mar 2025).

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 $Q1=0$6. Because object contributions are local—$Q1=0$7 when $Q1=0$8—the posterior has a local Markov structure. The variational family uses a $Q1=0$9-color checkerboard partition with rank

$\mathrm{rank}(Q)=n-1$0

and factorizes as

$\mathrm{rank}(Q)=n-1$1

For a posterior locality radius $\mathrm{rank}(Q)=n-1$2 and tile size $\mathrm{rank}(Q)=n-1$3, the minimal number of colors is $\mathrm{rank}(Q)=n-1$4. Within each tile, an auxiliary permutation $\mathrm{rank}(Q)=n-1$5 gives

$\mathrm{rank}(Q)=n-1$6

This yields a variational graph that mirrors the posterior conditional independencies while preserving GPU parallelism within each color class (Regier, 3 Oct 2025).

Latent Memory Palace adopts a variable-length discrete autoregressive latent sequence rather than a tiled field. With observation $\mathrm{rank}(Q)=n-1$7, action chunk $\mathrm{rank}(Q)=n-1$8, latent vocabulary $\mathrm{rank}(Q)=n-1$9, and termination time $K$0, the model uses

$K$1

with decoder $K$2. 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 (Zhu et al., 9 Jul 2026).

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: $K$3

$K$4

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 (Thakur, 8 Apr 2026).

Astronomical SAVD is trained not by an ELBO on observed data but by Neural Posterior Estimation under the expected forward KL. The objective is

$K$5

with unbiased gradient

$K$6

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 (Regier, 3 Oct 2025).

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,

$K$7

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 $K$8, so extra latent steps are useful only when they improve the reconstruction mean enough to offset tighter variance (Zhu et al., 9 Jul 2026).

For variational autoregressive networks in statistical mechanics, the objective is variational free energy rather than a standard data ELBO. With target Boltzmann distribution $K$9, minimizing $D=\mathrm{diag}(d_1,\dots,d_n)$00 is equivalent to minimizing

$D=\mathrm{diag}(d_1,\dots,d_n)$01

The natural-gradient update uses the Fisher information

$D=\mathrm{diag}(d_1,\dots,d_n)$02

and, in low-rank sample space, becomes

$D=\mathrm{diag}(d_1,\dots,d_n)$03

This changes the dominant inversion from parameter dimension $D=\mathrm{diag}(d_1,\dots,d_n)$04 to batch size $D=\mathrm{diag}(d_1,\dots,d_n)$05, giving complexity cubic in batch size rather than in model parameters (Liu et al., 2024).

Missing-data spatial models often require hybridization. In SEMs under MAR or MNAR, HVB factorizes as $D=\mathrm{diag}(d_1,\dots,d_n)$06 or its MNAR analogue and draws $D=\mathrm{diag}(d_1,\dots,d_n)$07 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 $D=\mathrm{diag}(d_1,\dots,d_n)$08-SEM variants (Wijayawardhana et al., 2024, Wijayawardhana et al., 29 May 2025). 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 (Gefang et al., 2022).

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 $D=\mathrm{diag}(d_1,\dots,d_n)$09 and MAE $D=\mathrm{diag}(d_1,\dots,d_n)$10, compared with Exact REML $D=\mathrm{diag}(d_1,\dots,d_n)$11, $D=\mathrm{diag}(d_1,\dots,d_n)$12, MLE $D=\mathrm{diag}(d_1,\dots,d_n)$13, $D=\mathrm{diag}(d_1,\dots,d_n)$14, and INLA $D=\mathrm{diag}(d_1,\dots,d_n)$15, $D=\mathrm{diag}(d_1,\dots,d_n)$16. 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 $D=\mathrm{diag}(d_1,\dots,d_n)$17 (Thakur, 8 Apr 2026).

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 $D=\mathrm{diag}(d_1,\dots,d_n)$18 for $D=\mathrm{diag}(d_1,\dots,d_n)$19 to $D=\mathrm{diag}(d_1,\dots,d_n)$20 for $D=\mathrm{diag}(d_1,\dots,d_n)$21, and F1 improves from approximately $D=\mathrm{diag}(d_1,\dots,d_n)$22 to $D=\mathrm{diag}(d_1,\dots,d_n)$23. In the crowded M2 setting, log-likelihood improves from $D=\mathrm{diag}(d_1,\dots,d_n)$24 to $D=\mathrm{diag}(d_1,\dots,d_n)$25, and F1 improves from approximately $D=\mathrm{diag}(d_1,\dots,d_n)$26 to $D=\mathrm{diag}(d_1,\dots,d_n)$27. The paper repeatedly attributes these gains to reduced border artifacts and better calibration near tile boundaries (Regier, 3 Oct 2025).

In control, LMP-$D=\mathrm{diag}(d_1,\dots,d_n)$28 and LMP-$D=\mathrm{diag}(d_1,\dots,d_n)$29 connect spatial retrieval to adaptive test-time compute. On DROID zero-shot tasks, LMP-$D=\mathrm{diag}(d_1,\dots,d_n)$30 reports block-bowl $D=\mathrm{diag}(d_1,\dots,d_n)$31 versus $D=\mathrm{diag}(d_1,\dots,d_n)$32 for DP and marker-mug $D=\mathrm{diag}(d_1,\dots,d_n)$33 versus $D=\mathrm{diag}(d_1,\dots,d_n)$34. On finetuned DROID tasks, it reports peg-hole $D=\mathrm{diag}(d_1,\dots,d_n)$35 versus $D=\mathrm{diag}(d_1,\dots,d_n)$36, and clean-table $D=\mathrm{diag}(d_1,\dots,d_n)$37, $D=\mathrm{diag}(d_1,\dots,d_n)$38, $D=\mathrm{diag}(d_1,\dots,d_n)$39 values of $D=\mathrm{diag}(d_1,\dots,d_n)$40 versus $D=\mathrm{diag}(d_1,\dots,d_n)$41, $D=\mathrm{diag}(d_1,\dots,d_n)$42 versus $D=\mathrm{diag}(d_1,\dots,d_n)$43, and $D=\mathrm{diag}(d_1,\dots,d_n)$44 versus $D=\mathrm{diag}(d_1,\dots,d_n)$45. On LIBERO-90 multitask, average performance is $D=\mathrm{diag}(d_1,\dots,d_n)$46 versus $D=\mathrm{diag}(d_1,\dots,d_n)$47, and bottom-10 tasks are $D=\mathrm{diag}(d_1,\dots,d_n)$48 versus $D=\mathrm{diag}(d_1,\dots,d_n)$49. The paper further reports a negative correlation $D=\mathrm{diag}(d_1,\dots,d_n)$50 between average latent steps and KNN action variance, and interprets this as linking adaptive compute to irreducible uncertainty (Zhu et al., 9 Jul 2026).

Autoencoded spatial priors in Fay–Herriot models are motivated primarily by scalability. In Missouri simulations, SMS-FH requires about $D=\mathrm{diag}(d_1,\dots,d_n)$51 hours per simulation and VSMS-FH about $D=\mathrm{diag}(d_1,\dots,d_n)$52 minutes, while GMS-FH requires about $D=\mathrm{diag}(d_1,\dots,d_n)$53 hours and VGMS-FH about $D=\mathrm{diag}(d_1,\dots,d_n)$54 minutes. In California, the paper states that GMS-FH would require more than $D=\mathrm{diag}(d_1,\dots,d_n)$55 days for one simulation, whereas VGMS-FH requires about $D=\mathrm{diag}(d_1,\dots,d_n)$56 hours and VSMS-FH about $D=\mathrm{diag}(d_1,\dots,d_n)$57 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 (Wang et al., 18 Mar 2025).

For SEMs with missing data, the main practical lesson is that variational structure alone is not always sufficient. Under MAR with $D=\mathrm{diag}(d_1,\dots,d_n)$58 and $D=\mathrm{diag}(d_1,\dots,d_n)$59 missing, JVB gives $D=\mathrm{diag}(d_1,\dots,d_n)$60 and $D=\mathrm{diag}(d_1,\dots,d_n)$61, whereas HVB-G gives $D=\mathrm{diag}(d_1,\dots,d_n)$62 and $D=\mathrm{diag}(d_1,\dots,d_n)$63. Under MNAR with the same dimensions and missingness, HVB-AllB and HVB-3B give $D=\mathrm{diag}(d_1,\dots,d_n)$64 and $D=\mathrm{diag}(d_1,\dots,d_n)$65, versus JVB $D=\mathrm{diag}(d_1,\dots,d_n)$66. 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 (Wijayawardhana et al., 2024, Wijayawardhana et al., 29 May 2025).

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 (Zhu et al., 9 Jul 2026, Regier, 3 Oct 2025, Thakur, 8 Apr 2026).

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 $D=\mathrm{diag}(d_1,\dots,d_n)$67 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 (Zhu et al., 9 Jul 2026, Regier, 3 Oct 2025, Wijayawardhana et al., 2024, Gefang et al., 2022).