Multivariate Autoregressive Training Paradigm

• The multivariate autoregressive training paradigm is a framework that develops strategies exploiting autoregressive structures to capture temporal and cross-dimensional dependencies.
• It decomposes dynamics into linear autoregressive and nonlinear cross-relation components, enabling specialized optimization and error correction across forecasting horizons.
• Empirical studies demonstrate reduced forecast errors and enhanced stability across linear, nonlinear, and probabilistic settings in diverse application domains.

Searching arXiv for the key papers to ground the article in current arXiv records. The multivariate autoregressive training paradigm denotes a class of training strategies for models that forecast a multivariate process by making autoregressive structure explicit at the level of parameterization, optimization, or rollout. In recent work, this has been formulated either as repeated application of a single-step multivariate predictor over future horizons, as in MVAR, or as an explicit decomposition of a multivariate transition operator into autoregressive and cross-relation components, as in AltTS (Fan et al., 16 Jul 2025, Yuan et al., 12 Feb 2026). Across the literature, the paradigm appears in linear, nonlinear, probabilistic, residual, directional, and manifold-valued settings, but its unifying feature is that the training procedure is designed around the fact that multivariate prediction errors propagate through time and across dimensions rather than remaining isolated within independent channels (Goel et al., 2017, Rasul et al., 2021).

1. Formal definitions and operator viewpoint

A common formalization starts from multivariate long-term time series forecasting. In AltTS, the lookback window is $X_t \in \mathbb{R}^{D\times L}$, the target horizon is $Y_{t+1} \in \mathbb{R}^{D\times H}$, and the model is expressed through a transition matrix $\mathcal{F} = (f_{ij})_{D\times D}$ with blocks $f_{ij}:\mathbb{R}^L \to \mathbb{R}^H$. Using the paper’s “apply-then-sum” operator,

$$\mathcal{F} * X_t := \Big(\sum_{j=1}^D f_{1j}(X_j),\dots,\sum_{j=1}^D f_{Dj}(X_j)\Big),$$

the ideal transition equation is

$\mathcal{F} * X_t + \varepsilon_{t+1} = Y_{t+1},$

and training minimizes

$$\mathcal{L}(\theta)=\sum_t \ell(\hat{\mathcal{F}}_\theta * X_t,\; Y_{t+1}).$$

This formulation isolates diagonal autoregressive effects $f_{ii}$ from off-diagonal cross-dimension effects $f_{ij}$, and it makes explicit that multivariate forecasting is not merely a collection of independent univariate problems (Yuan et al., 12 Feb 2026).

Other formulations preserve the same core logic while changing the state space. In MVAR, the multivariate state at time $t$ is $Y_{t+1} \in \mathbb{R}^{D\times H}$0, with $Y_{t+1} \in \mathbb{R}^{D\times H}$1 cities and $Y_{t+1} \in \mathbb{R}^{D\times H}$2 pollutants, and the multi-step mapping is defined by repeated application of a single-step model $Y_{t+1} \in \mathbb{R}^{D\times H}$3: $Y_{t+1} \in \mathbb{R}^{D\times H}$4 Without meteorology, the recursion is

$Y_{t+1} \in \mathbb{R}^{D\times H}$5

and with meteorology it is

$Y_{t+1} \in \mathbb{R}^{D\times H}$6

so the training target is inherently autoregressive and multivariate at every rollout step (Fan et al., 16 Jul 2025).

Probabilistic variants factorize the forecast distribution over time rather than directly over coordinates. TimeGrad models the conditional future law as

$Y_{t+1} \in \mathbb{R}^{D\times H}$7

where each factor is a multivariate distribution over $Y_{t+1} \in \mathbb{R}^{D\times H}$8. The model is thus autoregressive in time but not factorized across variables within a time slice (Rasul et al., 2021). This suggests that “autoregressive” in the multivariate literature refers less to a single architectural template than to a family of factorizations and rollout regimes that preserve temporal causality while modeling dependence across dimensions.

2. Structural decompositions of multivariate dependence

A central theme is decomposition: separate the easy, stable, within-series dynamics from the harder, higher-variance, cross-series component. AltTS states this directly by splitting the transition operator as

$Y_{t+1} \in \mathbb{R}^{D\times H}$9

with $\mathcal{F} = (f_{ij})_{D\times D}$0 and $\mathcal{F} = (f_{ij})_{D\times D}$1. The AR path is a channel-independent linear forecaster with RevIN normalization and an $\mathcal{F} = (f_{ij})_{D\times D}$2 regularizer, while the CR path is a Transformer over variable tokens using Cross-Relation Self-Attention, which masks the diagonal so that a variable cannot attend to itself. The model therefore enforces role specialization structurally: autoregression remains explicit and linear, while cross-dimension interactions are nonlinear and masked to remain genuinely cross-variable (Yuan et al., 12 Feb 2026).

Earlier hybrid models instantiated related decompositions without the same optimization argument. R2N2 first fits a simple linear multivariate autoregressive model, typically VAR-1, then trains an LSTM on the residual multivariate time series and predicts

$\mathcal{F} = (f_{ij})_{D\times D}$3

The linear model captures short-lag linear dependencies and the recurrent model captures residual nonlinear structure and longer-term effects (Goel et al., 2017). The autoregressive convolutional recurrent neural network of 2019 similarly combines a multi-scale convolutional and GRU-based nonlinear path with a linear autoregressive shortcut

$\mathcal{F} = (f_{ij})_{D\times D}$4

where the linear term uses a regression input window of 5 previous steps and shares weights and biases across variables (Maggiolo et al., 2019). In both cases, the final prediction is additive, and the linear autoregressive component functions as a stable baseline while the nonlinear module supplies corrective structure.

Recent high-dimensional work retains the same separation principle but introduces stronger parameter sharing. RRNAR represents the transition as

$\mathcal{F} = (f_{ij})_{D\times D}$5

with

$\mathcal{F} = (f_{ij})_{D\times D}$6

The network side is reduced to two scalars, while the variable side is low rank; this couples known topology with a learnable low-rank variable subspace and yields a parameter count of $\mathcal{F} = (f_{ij})_{D\times D}$7 rather than $\mathcal{F} = (f_{ij})_{D\times D}$8 (Lyu et al., 4 Jan 2026). A plausible implication is that multivariate autoregressive training increasingly treats parameterization as an optimization device: structural decomposition is used not only for interpretability but also to control variance and dimensionality.

3. Optimization schedules as part of the model

One of the most explicit claims in this literature is that training schedule is itself part of the model class. AltTS argues that joint optimization of AR and CR parameters produces gradient entanglement because the aggregate residual

$\mathcal{F} = (f_{ij})_{D\times D}$9

mixes diagonal and off-diagonal errors, so the AR gradient becomes

$f_{ij}:\mathbb{R}^L \to \mathbb{R}^H$0

The cross-relation term contaminates the autoregressive update, increases covariance, and destabilizes long-horizon learning. The proposed response is block alternating optimization with two independent AMSGrad optimizers: first update $f_{ij}:\mathbb{R}^L \to \mathbb{R}^H$1 with $f_{ij}:\mathbb{R}^L \to \mathbb{R}^H$2 fixed, then update $f_{ij}:\mathbb{R}^L \to \mathbb{R}^H$3 with $f_{ij}:\mathbb{R}^L \to \mathbb{R}^H$4 fixed. The paper reports that AR is always updated before CR and notes a fixed rhythm such as 10 AR updates followed by 2 CR updates per mini-batch (Yuan et al., 12 Feb 2026).

MVAR treats training-time rollout consistency as the decisive issue. Its Multivariate Autoregressive Training Paradigm trains a single-step operator by free-running autoregression rather than teacher forcing. The model takes only 2 historical time steps, forecasts up to 20 steps, and uses the Step Weighted loss

$f_{ij}:\mathbb{R}^L \to \mathbb{R}^H$5

where $f_{ij}:\mathbb{R}^L \to \mathbb{R}^H$6 is a monotonically decreasing, equally spaced sequence with maximum value 5 and minimum 0.1. The paper explicitly states that teacher forcing is not used; during training, predictions are always fed back as inputs to the next step (Fan et al., 16 Jul 2025).

BPTT-SA addresses the same mismatch in recurrent forecasting. Standard teacher-forced BPTT computes gradients under the data distribution, whereas autoregressive inference feeds back model predictions. BPTT-SA samples Bernoulli switches $f_{ij}:\mathbb{R}^L \to \mathbb{R}^H$7 and interpolates between teacher forcing and autoregressive propagation; its hidden-state gradient contains products of either the teacher-forcing factor $f_{ij}:\mathbb{R}^L \to \mathbb{R}^H$8 or the autoregressive factor $f_{ij}:\mathbb{R}^L \to \mathbb{R}^H$9, yielding

$$\mathcal{F} * X_t := \Big(\sum_{j=1}^D f_{1j}(X_j),\dots,\sum_{j=1}^D f_{Dj}(X_j)\Big),$$0

Early in training $$\mathcal{F} * X_t := \Big(\sum_{j=1}^D f_{1j}(X_j),\dots,\sum_{j=1}^D f_{Dj}(X_j)\Big),$$1, later $$\mathcal{F} * X_t := \Big(\sum_{j=1}^D f_{1j}(X_j),\dots,\sum_{j=1}^D f_{Dj}(X_j)\Big),$$2, and validation is computed with $$\mathcal{F} * X_t := \Big(\sum_{j=1}^D f_{1j}(X_j),\dots,\sum_{j=1}^D f_{Dj}(X_j)\Big),$$3, namely pure autoregressive inference (Vlachas et al., 2023).

Not all paradigms are simultaneous or end-to-end. R2N2 is explicitly sequential: fit VAR, freeze it, compute residuals, and then train the RNN on residual-augmented sequences (Goel et al., 2017). By contrast, the 2019 convolutional recurrent model trains its convolutional, recurrent, and linear AR components jointly under MSE (Maggiolo et al., 2019). Taken together, these works indicate that “training paradigm” in this area includes at least four distinct regimes: joint additive training, residual-stage training, alternating block optimization, and closed-loop free-running rollout training.

4. Probabilistic, bidirectional, and non-Euclidean extensions

The paradigm generalizes beyond deterministic Euclidean regression. TimeGrad preserves the temporal autoregressive factorization but replaces a simple emission head with a conditional diffusion model. At each step $$\mathcal{F} * X_t := \Big(\sum_{j=1}^D f_{1j}(X_j),\dots,\sum_{j=1}^D f_{Dj}(X_j)\Big),$$4, the RNN hidden state $$\mathcal{F} * X_t := \Big(\sum_{j=1}^D f_{1j}(X_j),\dots,\sum_{j=1}^D f_{Dj}(X_j)\Big),$$5 conditions a diffusion loss

$$\mathcal{F} * X_t := \Big(\sum_{j=1}^D f_{1j}(X_j),\dots,\sum_{j=1}^D f_{Dj}(X_j)\Big),$$6

so the model learns one-step multivariate predictive distributions and then rolls them forward autoregressively at inference (Rasul et al., 2021). This suggests that autoregressive training need not imply a simple Gaussian head; it can be combined with expressive score-based generative modeling.

COPAR replaces linear dependence with a vine-copula factorization. In the bivariate case, serial dependence is represented by $$\mathcal{F} * X_t := \Big(\sum_{j=1}^D f_{1j}(X_j),\dots,\sum_{j=1}^D f_{Dj}(X_j)\Big),$$7, conditional serial dependence in the second series by $$\mathcal{F} * X_t := \Big(\sum_{j=1}^D f_{1j}(X_j),\dots,\sum_{j=1}^D f_{Dj}(X_j)\Big),$$8, and cross-series dependence by $$\mathcal{F} * X_t := \Big(\sum_{j=1}^D f_{1j}(X_j),\dots,\sum_{j=1}^D f_{Dj}(X_j)\Big),$$9 and $\mathcal{F} * X_t + \varepsilon_{t+1} = Y_{t+1},$0. A COPAR($\mathcal{F} * X_t + \varepsilon_{t+1} = Y_{t+1},$1) model sets all copulas with lag greater than $\mathcal{F} * X_t + \varepsilon_{t+1} = Y_{t+1},$2 to the independence copula and therefore uses $\mathcal{F} * X_t + \varepsilon_{t+1} = Y_{t+1},$3 non-independence copulas in the bivariate case (Brechmann et al., 2012). The autoregressive paradigm here is expressed through lag truncation and sequential pair-copula estimation rather than through neural rollout.

MIM extends autoregression along direction rather than variable axis. It jointly trains left-to-right and right-to-left factorizations,

$\mathcal{F} * X_t + \varepsilon_{t+1} = Y_{t+1},$4

with a token-wise total variation alignment term and a bidirectional inference procedure that lets both directional models “meet in the middle” (Nguyen et al., 2023). Although this is not a multivariate time-series paper, it broadens the meaning of autoregressive training by showing that multiple conditional factorizations of the same joint distribution can be trained jointly and aligned.

The same pattern appears in non-Euclidean time series. The Wasserstein multivariate autoregressive model centers distribution-valued observations by Fréchet means and works in the tangent space at the Lebesgue measure. Its centered quantile dynamics take the form

$\mathcal{F} * X_t + \varepsilon_{t+1} = Y_{t+1},$5

with simplex constraints

$\mathcal{F} * X_t + \varepsilon_{t+1} = Y_{t+1},$6

to ensure valid quantile maps (Jiang et al., 2022). This is still multivariate autoregression, but training occurs in a Wasserstein tangent geometry and the estimator is a constrained quadratic program rather than a standard least-squares fit.

5. High-dimensional, structured, and statistical formulations

A prominent branch of the paradigm is devoted to dimensionality control in high-dimensional settings. TVP-MAI reduces a large TVP-VAR-SV into a small number of indices $\mathcal{F} * X_t + \varepsilon_{t+1} = Y_{t+1},$7, then models

$\mathcal{F} * X_t + \varepsilon_{t+1} = Y_{t+1},$8

with state-space recursion

$\mathcal{F} * X_t + \varepsilon_{t+1} = Y_{t+1},$9

Training mixes a switching algorithm for $$\mathcal{L}(\theta)=\sum_t \ell(\hat{\mathcal{F}}_\theta * X_t,\; Y_{t+1}).$$0, forgetting-factor Kalman filtering for $$\mathcal{L}(\theta)=\sum_t \ell(\hat{\mathcal{F}}_\theta * X_t,\; Y_{t+1}).$$1, EWMA volatility updates, and Dynamic Model Selection or Dynamic Model Averaging over $$\mathcal{L}(\theta)=\sum_t \ell(\hat{\mathcal{F}}_\theta * X_t,\; Y_{t+1}).$$2 choices (Cubadda et al., 2022). The paper emphasizes that this avoids the computational burden of large MCMC-based TVP-VAR-SV models.

TVART assumes a windowed VAR(1),

$$\mathcal{L}(\theta)=\sum_t \ell(\hat{\mathcal{F}}_\theta * X_t,\; Y_{t+1}).$$3

stacks $$\mathcal{L}(\theta)=\sum_t \ell(\hat{\mathcal{F}}_\theta * X_t,\; Y_{t+1}).$$4 into a tensor $$\mathcal{L}(\theta)=\sum_t \ell(\hat{\mathcal{F}}_\theta * X_t,\; Y_{t+1}).$$5, and imposes a CP decomposition

$$\mathcal{L}(\theta)=\sum_t \ell(\hat{\mathcal{F}}_\theta * X_t,\; Y_{t+1}).$$6

The regularized objective is

$$\mathcal{L}(\theta)=\sum_t \ell(\hat{\mathcal{F}}_\theta * X_t,\; Y_{t+1}).$$7

with $$\mathcal{L}(\theta)=\sum_t \ell(\hat{\mathcal{F}}_\theta * X_t,\; Y_{t+1}).$$8 equal to temporal total variation or spline smoothing. Training uses alternating minimization over the three factor matrices, and the paper proves that any stationary point of the algorithm is a local minimum (Harris et al., 2019).

In multivariate ordinal regression with repeated measurements, the autoregressive structure is moved to latent errors: $$\mathcal{L}(\theta)=\sum_t \ell(\hat{\mathcal{F}}_\theta * X_t,\; Y_{t+1}).$$9 The full $f_{ii}$0 covariance $f_{ii}$1 is induced by $f_{ii}$2 and $f_{ii}$3, but estimation avoids $f_{ii}$4-dimensional integrals by maximizing a pairwise composite likelihood over bivariate ordinal probabilities (Vana-Gür, 2024). Similarly, the biomedical signal-processing framework introduces an overparameterized loss

$f_{ii}$5

then alternates between parameter estimation and state reconstruction via explicit linear solves, including multivariate AR(1) EEG connectivity with $f_{ii}$6 channels (Haderlein et al., 2023). These statistical formulations show that the training paradigm can be expressed as composite likelihood or alternating latent-state estimation just as naturally as it can be expressed through deep learning optimizers.

6. Empirical behavior, strengths, and limitations

Across the cited papers, the main empirical claim is not that one universal architecture dominates, but that training design materially affects multivariate forecasting quality. AltTS reports that across 56 dataset-horizon-metric combinations it is best or second best in 49 cases, with particularly strong gains on ETTh1/2, ETTm1/2, and Weather at long horizons; in the AO ablation, ETTh1 at $f_{ii}$7 improves from $f_{ii}$8 without AO to $f_{ii}$9 with AO (Yuan et al., 12 Feb 2026). MVAR reports that a model trained with only 2 past steps can outperform state-of-the-art baselines that use 20-step inputs, and that meteorology-coupled variants further improve PM$f_{ij}$0 RMSE, for example from 29.66 to 28.11 at 1–24h and from 34.10 to 31.29 at 97–120h when ERA5 is used (Fan et al., 16 Jul 2025).

The same pattern appears outside deterministic forecasting. TimeGrad reports state-of-the-art multivariate probabilistic forecasting on datasets with thousands of correlated dimensions and attributes this to combining an autoregressive temporal encoder with a diffusion-based multivariate conditional head (Rasul et al., 2021). In offline RL, autoregressive dynamics models outperform diagonal-Gaussian baselines in held-out negative log-likelihood on every DM Control task and achieve a median OPE rank correlation of 0.90, a median normalized Regret@5 of 0.04, and a median absolute error of 9.71, surpassing model-based diagonal baselines and model-free OPE methods (Zhang et al., 2021). In MIM, bidirectional autoregressive pretraining reduces OpenWebText perplexity from 11.92 to 9.54 at 2.7B and improves HumanEval Infilling pass@1 from 22.8% to 26.3% relative to FIM at the same scale (Nguyen et al., 2023). These results suggest that the autoregressive training paradigm has become a general methodological principle rather than a domain-specific recipe.

The literature also identifies recurrent limitations. Linear AR backbones can be too restrictive when nonlinear within-series effects dominate (Yuan et al., 12 Feb 2026). Free-running autoregressive training can accumulate error on spike-driven targets such as PM$f_{ij}$1 (Fan et al., 16 Jul 2025). Diffusion-based autoregressive heads can be slow because each forecast step requires multiple reverse-process evaluations (Rasul et al., 2021). Autoregressive dynamics models in RL require $f_{ij}$2 forward passes across state dimensions at inference (Zhang et al., 2021). Composite-likelihood and sandwich-variance procedures can underestimate uncertainty for thresholds and may need bootstrap or jackknife corrections (Vana-Gür, 2024). These are not contradictions of the paradigm; they indicate that explicit autoregression does not remove the classical trade-offs among bias, variance, computational cost, and model misspecification.

A persistent misconception is that multivariate autoregressive training is simply “teacher forcing on a multivariate output.” The surveyed papers do not support that reduction. Some methods are explicitly free-running and avoid teacher forcing during training (Fan et al., 16 Jul 2025); some alternate among blocks to suppress gradient contamination (Yuan et al., 12 Feb 2026); some decompose residuals sequentially (Goel et al., 2017); some estimate latent states and parameters jointly (Haderlein et al., 2023); and some factorize the conditional law through copulas, diffusion chains, or bidirectional sequence factorizations (Brechmann et al., 2012, Rasul et al., 2021, Nguyen et al., 2023). A more accurate characterization is that the paradigm treats autoregressive dependence as a modeling object and the training schedule as a statistical object of equal importance.