FAMoE-ST: Hierarchically Frozen Attention with Hybrid-Memory Experts for Traffic Flow Forecasting

1. Introduction

Traffic flow forecasting estimates future traffic states from historical roadside-sensor observations and supports signal control, route planning, congestion warning, and traffic guidance. The task is difficult because traffic dynamics are coupled across both space and time. Dependencies between nearby and distant sensors change as congestion propagates, while daily and weekly periodicity interacts with incidents, weather, and localized demand surges. A reliable multi-node, multi-step forecaster must therefore represent dynamic spatial interactions, multiscale temporal variation, and heterogeneous node behavior, and must remain stable at long horizons and under high-flow conditions.

Recent task-specific models have improved traffic propagation modeling through graph neural differential equations, spatiotemporal fusion, and dynamic graph decoupling [1,2,3,4]. STID, STEP, PDFormer, and STAEformer further strengthen multi-step forecasting through identity embeddings, self-supervised pretraining, propagation-delay modeling, and adaptive spatiotemporal embeddings [5,6,7,8], whereas TESTAM uses a mixture of experts to accommodate changes in traffic regimes [9]. These methods provide effective inductive biases, but many rely on task-specific graph construction or specialized architectural components, which can complicate adaptation across datasets and operating conditions.

A complementary line of research repurposes large pretrained or generally structured sequence models for traffic and time-series forecasting. ST-LLM introduces a GPT-style backbone into traffic prediction [10], while One Fits All, Time-LLM, zero-shot language-model forecasting, LLM4TS, TimesFM, Chronos, and MOMENT explore input mapping, reprogramming, or foundation-model pretraining [11,12,13,14,15,16,17]. PatchTST and TimesNet likewise demonstrate the value of generic sequence architectures for long-range and multiscale patterns [18,19]. However, language tokens and traffic measurements differ in semantics, distribution, and dependency structure. In this study, GPT-2 supplies the initial Transformer weights and block structure; beneficial transfer of linguistic knowledge is not assumed. Establishing such transfer would require a randomly initialized isomorphic backbone under the same training protocol. Full adaptation provides flexibility but increases optimization cost and the risk of unstable updates, whereas complete freezing can underfit the target domain. Parameter-efficient methods such as LoRA and P-Tuning address a related trade-off by restricting the trainable parameter space [20,21]. In parallel, conditional-computation models expand representational capacity through sparse or soft expert assignment [22,23,24,25], yet routers driven only by the current input can be sensitive to short-lived noise and imbalanced token allocation. A unified design that controls backbone adaptation and combines instantaneous traffic states with cross-batch routing references remains underexplored.

To address these limitations, we propose FAMoE-ST, which combines depth-controlled adaptation of a GPT-style spatial Transformer with hybrid-memory expert routing. The principal contributions are as follows:

• A hierarchical freezing strategy fixes the lower attention and feed-forward blocks while adapting the upper blocks, thereby limiting optimization to the depth at which traffic-specific spatial representations and conditional transformations are learned.
• An ordering-neutral spatial attention design treats sensors as tokens, replaces GPT-2's causal mask with bidirectional all-to-all attention, and assigns the shared zero GPT position ID to every token so that the available spatial context does not depend on arbitrary node numbering.
• A hybrid-memory mixture-of-experts module combines an input-responsive linear router with a learnable cross-batch memory router and Top-2 sparse gating. Experiments on PEMS04 and PEMS08, including baseline comparisons, component ablations, data-ratio sensitivity tests, routing diagnostics, and a parameter-matched dense control, quantify the complementary effects of freezing, expert capacity, and memory-assisted routing while separating routing gains from raw parameter growth.

2. Materials and Methods

2.1. Problem Formulation and Overall Architecture

Consider a traffic network with $N$ sensor nodes. For sample $b$, $X$ denotes the historical observation tensor over the $P$ time steps preceding forecast origin t and has dimensions $P\times N\times C_x$; $Y$ denotes the ground-truth traffic states over the subsequent $S$ time steps and has dimensions $S\times N\times C_y$. Here, $C_x$ includes the traffic variable and available temporal covariates, whereas $C_y$is the number of target variables. The objective is to learn a parameterized mapping that directly transforms the historical window into a multi-step forecast:

$${\widehat{Y}}_{(t+1:t+S)}={\mathcal{F}}_{\mathsf{\Theta }}(X_{(t-P+1:t)}),\widehat{Y}\in {\mathbb{R}}^{S\times N\times C_y}$$

(1)

where X is the historical observation tensor,$\widehat{Y}$ is the predicted traffic-state tensor, t is the forecast origin, $P$ and $S$ are the input and output window lengths, $N$ is the number of sensor nodes, $C_y$ is the target dimension, $\mathsf{\Theta }$ is the set of trainable parameters, ${\mathcal{F}}_{\mathsf{\Theta }}$ is the parameterized forecasting function, and ℝ is the real-number domain.

In the 12-step setting considered here, $P=S=12$ and $C_y$ = 1. Each sensor is represented by one token that summarizes its complete P-step history and associated covariates. The N sensor tokens form the sequence processed by the backbone; thus, the Transformer's sequence axis is spatial rather than temporal. This construction allows cross-node interaction and node-specific expert assignment within a common representation space.

This formulation presents two central challenges. First, inter-node dependencies vary with traffic conditions, so fixed spatial relations cannot fully describe dynamic propagation. Second, sensors and time periods exhibit distinct operating regimes that a single shared feed-forward network may represent inefficiently. FAMoE-ST addresses these challenges with a spatiotemporal node tokenizer, ordering-neutral bidirectional spatial attention, hierarchically restricted adaptation, and hybrid-memory experts.

As illustrated in Figure 1, the model first converts each sensor's historical sequence and context into a node token. A Transformer backbone initialized with GPT-2 weights then models global cross-node dependencies. Lower blocks remain fixed, upper-layer attention is trainable, and the upper feed-forward networks are replaced by hybrid-memory expert modules that select token-specific transformations. A lightweight prediction head produces all forecast horizons in one pass. The end-to-end mapping is

$$H^0={\mathsf{\Phi }}_{tok}(X),H^L={\mathsf{\Phi }}_{FA-HMoE}(H^0),\widehat{Y}={\mathsf{\Phi }}_{pred}(H^L)$$

(2)

where $H^0$and $H^L$ denote the backbone input and the node representations after L layers, respectively; ${\mathsf{\Phi }}_{tok}$ is the spatiotemporal node tokenizer; ${\mathsf{\Phi }}_{FA-HMoE}$ is the hierarchically frozen-attention backbone with hybrid-memory experts; ${\mathsf{\Phi }}_{pred}$is the prediction head; and L is the number of backbone layers. $X$ and $\widehat{Y}$ retain the definitions given for Equation (1).

The tokenizer ${\mathsf{\Phi }}_{tok}$ maps continuous traffic histories and contextual variables into the hidden space expected by the Transformer. The backbone ${\mathsf{\Phi }}_{FA-HMoE}$ fixes the attention and feed-forward parameters in lower layers while jointly updating attention, routing, and expert parameters in upper layers. The head ${\mathsf{\Phi }}_{pred}$ projects the final node representations onto $S$ forecast steps. This depth-wise organization controls optimization cost and adaptation capacity, while conditional computation accommodates heterogeneity across nodes and traffic regimes.

At the input, the spatiotemporal node tokenizer projects the historical window, periodic context, and node identity into a shared hidden space. For node n, the model constructs traffic-state, time-of-day/day-of-week, and node-identity embeddings separately. Their joint encoding is

$$e_{b,n}^x=W_{x}vec(X_{b,:,n,:})+b_x,e_b^t=E_{day}\left[{\tau }_b\right]+E_{week}\left[{\omega }_b\right],e_n^s=E_{node}\left[n\right]$$

$${\mathrm{h}}_{b,n}^0=W_f[e_{b,n}^x\Vert e_b^t\Vert e_n^s]+b_f+p_n,H_b^0=[{\mathrm{h}}_{b,1}^0,\dots ,{\mathrm{h}}_{b,N}^0]\in {\mathbb{R}}^{N\times d_g}$$

(3)

where $e^x$, $e^t$, and $e^s$ are the traffic-state, temporal-position, and node-identity embeddings for sample $b$ and node $n$; $W_x$and $b_x$ are the traffic projection parameters; vec denotes vectorization; $E_{day}$, $E_{week}$, and $E_{node}$ are the time-of-day, day-of-week, and node embedding tables; ${\tau }_b$ and ${\omega }_b$ are the corresponding temporal indices; $W_f$ and $b_f$ are the fusion projection parameters; pn is the GPT positional term; h⁰ is the initial node representation; $H^0$stacks all $N$ node tokens in a sample; $d_g$ is the Transformer hidden dimension; and $\Vert$ denotes feature concatenation. In the implementation, $p_n=p_0$ for every node because all node tokens use the shared zero GPT position ID. Sensor identity is therefore carried by Enode rather than by an order-dependent GPT position.

The traffic projection is shared across nodes, the periodic embeddings identify positions within the day and week, and$E_{node}$ provides a learned sensor identity. Concatenation occurs only along the feature dimension and introduces no cross-node information before the backbone. After fusion into the Transformer hidden space, every sensor token may attend to every other sensor token. Because the mask is all-visible and every token uses the same zero GPT position ID, the available spatial context is independent of the sensor's index in the NPZ file; Enode, rather than sequence position, preserves sensor identity.

2.2. Hierarchically Frozen Attention

Adapting all backbone layers is costly and can produce unnecessary parameter updates when the target datasets are comparatively small, whereas complete freezing can prevent the backbone from learning traffic-specific interactions. We therefore restrict optimization by depth. Given backbone depth L, the numbers U and F of trainable upper and frozen lower layers, respectively, together with the layer-wise trainability indicator, are defined as follows:

$$U=\mathrm{m}\mathrm{a}\mathrm{x}\{1,\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}({\displaystyle \frac{L}{3}})\},F=L-U$$

$$m_{\mathcal{l}}=\left\{\begin{array}{ccc}0,& 1\le \mathcal{l}\le F& \mathrm{f}\mathrm{r}\mathrm{o}\mathrm{z}\mathrm{e}\mathrm{n};\\ 1,& F<\mathcal{l}\le L& \mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{b}\mathrm{l}\mathrm{e}.\end{array}\right.$$

$${\mathcal{S}}_M=\{F+1,\dots ,L\}=\left\{\mathrm{5,6}\right\}(L=6)$$

(4)

where L is the total number of backbone layers; U and F are the numbers of trainable upper layers and frozen lower layers, respectively; round and max denote the rounding and maximum operations; mℓ is the trainability indicator for layer ℓ, with $m_{\mathcal{l}}$ = 0 for a frozen layer and $m_{\mathcal{l}}$ = 1 for a trainable layer; and SM is the set of layers equipped with hybrid-memory experts. For L = 6, the rule gives U = 2, F = 4, and ${\mathcal{S}}_M$ = {5, 6}.

The rule freezes approximately two-thirds of the backbone. For L = 6, it gives U = 2 and F = 4: layers L1–L4 remain frozen, whereas L5–L6 are trainable. This allocation constrains the optimization search space while retaining upper-layer capacity for traffic-specific spatial adaptation. It also scales the freezing depth automatically when L changes, avoiding a separate manually selected layer count for each backbone depth. The backbone follows the pre-LayerNorm GPT-2 block structure. For block ℓ, multi-head attention is computed as

$$Z^{\mathcal{l}}={LN}_{\mathcal{l},1}(H^{\mathcal{l}-1})$$

$$Q_a^{\mathcal{l}}=Z^{\mathcal{l}}{W}_{Q,a}^{\mathcal{l}},K_a^{\mathcal{l}}=Z^{\mathcal{l}}{W}_{K,a}^{\mathcal{l}},V_a^{\mathcal{l}}=Z^{\mathcal{l}}{W}_{V,a}^{\mathcal{l}}$$

$$A_a^{\mathcal{l}}=softmax({\displaystyle \frac{Q_a^{\mathcal{l}}{K_a^{\mathcal{l}}}^T}{\sqrt{d_{\mathrm{h}}}}}+\Gamma )V_a^{\mathcal{l}}$$

$${MHA}_{\mathcal{l}}(Z^{\mathcal{l}})=Concat(A_1^{\mathcal{l}},\dots ,A_A^{\mathcal{l}})W_O^{\mathcal{l}}$$

(5)

where $Z^{\mathcal{l}}$ is the input after the first layer-normalization operation in layer ℓ; $LN$ denotes layer normalization; $Q$, $K$, and $V$ are the query, key, and value matrices; $W_Q$, $W_K$, $W_V$, and $W_O$ are the corresponding projections; a indexes an attention head and $A$ is the total number of heads; dh is the head dimension; Γ is the additive attention-bias matrix; $A_a^{\mathcal{l}}$ is the output of head a in layer ℓ; the superscript T denotes matrix transpose; and Concat and MHA denote head concatenation and multi-head attention. Sequence positions i and j index sensor tokens rather than forecast times. We set ${\Gamma }_{ij}$= 0 for every sensor pair (i, j), which yields bidirectional all-to-all spatial attention instead of the lower-triangular causal mask used in autoregressive language modeling. This operation introduces no future-target leakage because every token is formed exclusively from the historical input window. Both the attention and feed-forward sublayers retain residual connections:

$${\tilde{H}}^{\mathcal{l}}=H^{\mathcal{l}-1}+{MHA}_{\mathcal{l}}({LN}_{\mathcal{l},1}(H^{\mathcal{l}-1}))$$

$$H^{\mathcal{l}}={\tilde{H}}^{\mathcal{l}}+{\mathcal{G}}_{\mathcal{l}}({LN}_{\mathcal{l},2}({\tilde{H}}^{\mathcal{l}}))$$

$${\mathcal{G}}_{\mathcal{l}}=\left\{\begin{array}{cc}{FFN}_{\mathcal{l}}^{frozen},& 1\le \mathcal{l}\le F;\\ {HMoE}_{\mathcal{l}},& \mathcal{l}\in {\mathcal{S}}_M.\end{array}\right.$$

(6)

where $H^{\mathcal{l}-1}$ and Hℓ are the input and output hidden representations of layer ℓ; ${\tilde{H}}^{\mathcal{l}}$ is the intermediate state after the attention residual update; ${LN}_{\mathcal{l},1}$ and ${LN}_{\mathcal{l},2}$ are the two pre-LayerNorm operations; ${\mathcal{G}}_{\mathcal{l}}$is the layer's feed-forward transformation; FFNfrozen,ℓ is the frozen original feed-forward network; ${HMoE}_{\mathcal{l}}$ is the hybrid-memory expert module; $F$ is the number of frozen lower layers; and ${\mathcal{S}}_M$is the set of upper expert layers.

For ℓ ≤ $F$, the attention and original $FFN$ parameters are fixed, while LayerNorm parameters remain trainable to recalibrate activation statistics in the traffic domain. For ℓ > $F$, attention, LayerNorm, routing, and expert parameters are optimized jointly. Thus, freezing is used as a depth-control and optimization mechanism; the present experiments do not attribute its effect specifically to preservation of linguistic knowledge.

2.3. Hybrid-Memory Mixture-of-Experts

Mixture-of-experts (MoE) models expand representational capacity through conditional routing [22,23,24,25], and TESTAM shows that expert assignment can accommodate time-varying traffic patterns [9]. Nevertheless, a single $FFN$ applies the same transformation to every token, while an instantaneous router driven only by the current hidden state can overreact to short-term perturbations. In each of the upper U layers, we therefore replace the original GPT-2 multilayer perceptron with a Hybrid-Memory Sparse Experts (HMoE) module. For a mini-batch, each layer receives $T=BN$ node tokens and contains E candidate experts.

The experts are initialized as copies of the original GPT-2 multilayer perceptron weights and then differentiate as sparse routing exposes them to different token subsets. This copying provides a common initialization and does not by itself establish beneficial language-knowledge transfer. The memory module stores no raw traffic samples. Instead, learnable keys and expert labels act as routing prototypes that accumulate cross-batch pattern references. In parallel, the instantaneous branch maps each current token directly to an expert-probability vector and responds to short-term state changes:

$$p_i^{\mathrm{lin}}=\mathrm{softmax}(W_r{\mathrm{h}}_i),p_i^{\mathrm{lin}}\text{∈}{\mathsf{\Delta }}^{E-1}$$

(7)

where $hᵢ$ is the hidden representation of token i; $W_r$is the linear projection of the instantaneous router; ${\mathrm{p}}_{\mathrm{i}}^{\mathrm{lin}}$ is the resulting expert-probability vector; $E$ is the number of experts; and ${∆}^{E-1}$ is the $E-dimensional$ probability simplex.

Memory retrieval computes similarities between a normalized query and the memory keys and retains only the Km slots most relevant to the current token:

$${\tilde{q}}_i={\displaystyle \frac{tan\mathrm{h}(W_q{\mathrm{h}}_i)}{max({‖tan\mathrm{h}(W_q{\mathrm{h}}_i)‖}_2,\epsilon )}},{\tilde{k}}_m={\displaystyle \frac{k_m}{max({‖k_m‖}_2,\epsilon )}}$$

$$s_{i,m}={{\tilde{q}}_i}^T{\tilde{k}}_m, {\mathcal{T}}_i^m=TopK(s_i,K_m)$$

$$\overline{s_{i,m}}=\left\{\begin{array}{cc}{s}_{i,m},& m\in {\mathcal{T}}_i^m;\\ -\mathrm{\infty },& m\notin {\mathcal{T}}_i^m.\end{array}\right.$$

(8)

where ${\tilde{q}}_i$ and ${\tilde{k}}_m$ are the normalized query for token i and normalized key of memory slot m, respectively; Wq is the query projection; hᵢ is the token representation; kₘ is the learnable unnormalized key of slot m; ${‖\bullet ‖}_2$ is the Euclidean norm; and ε > 0 prevents division by zero. The raw similarity is $s_{i,m}={{\tilde{q}}_i}^T{\tilde{k}}_m$; $k_m$ is the number of retained memory slots; ${\mathcal{T}}_i^m$ is the set returned by $TopK(s_i,K_m)$; and s̄ᵢ,ₘ is the masked similarity, which equals $\overline{s_{i,m}}$ for selected slots and −∞ otherwise.

Normalizing similarities over the retained slots and mapping them to expert preferences yields the memory-routing distribution:

$${\alpha }_{i,m}={softmax}_m({\bar{s}}_{i,m}),{\rho }_m=softmax(r_m),p_i^{mem}={\sum }_{m=1}^M{\alpha }_{i,m}{\rho }_m$$

(9)

where$M$is the total number of memory slots and $m\in \{1,\dots ,M\}$ indexes a slot; ${\alpha }_{i,m}$ is the softmax-normalized weight of slot m for token i over the retained slots; rₘ is the learnable expert-logit vector attached to slot m; ${\rho }_m$ is its expert-preference distribution; and $p_i^{mem}$ is the memory-routing distribution.

Equation (9) defines a valid probability distribution. The memory branch can therefore supply a cross-batch routing reference by associating similar hidden patterns with expert preferences. The instantaneous and memory distributions are then fused by a convex combination with memory weight $\beta$:

$${\tilde{p}}_i=(1-\beta )p_i^{lin}+\beta p_i^{mem},\beta \in \left[\mathrm{0,1}\right]$$

$$p_i={\displaystyle \frac{{\tilde{p}}_i}{max({\displaystyle {\sum }_{e=1}^E}{\tilde{p}}_{i,e},\epsilon )}}$$

(10)

where ${\tilde{p}}_i$is the convex fusion of the instantaneous distribution $p_i^{lin}$ and memory distribution $p_i^{mem}$; $\beta \in [0,1]$is the memory-branch weight and $1-\beta$is the instantaneous-branch weight; $p_i$ is the numerically normalized fused expert distribution; $E$ is the number of experts; e indexes an expert; and ε prevents division by zero.

Because both branches are normalized, their convex combination is theoretically a valid probability distribution; the second normalization in Equation (10) compensates only for finite-precision error. The $K_e$ experts with the largest fused probabilities are selected, and their sparse gating weights are

$${\mathcal{T}}_i=TopK(p_i,K_e)$$

$$g_{i,e}=\left\{\begin{array}{cc}{\displaystyle \frac{p_{i,e}}{max(\sum _{j\in {\mathcal{T}}_i}{p}_{i,j},\epsilon )}},& e\in {\mathcal{T}}_i;\\ 0,& e\notin {\mathcal{T}}_i.\end{array}\right.$$

(11)

$$HMoE({\mathrm{h}}_i)={\sum }_{e=1}^E{g}_{i,e}{f}_e({\mathrm{h}}_i),{\mathcal{T}}_i=TopK(p_i,K_e)$$

(12)

where $K_e$ is the number of experts activated for each token; ${\mathcal{T}}_i$ is the selected expert set; $e$ and$j$index experts; $g_{i,e}$ is the gate weight assigned from token i to expert e after renormalization over Tᵢ; fₑ is the feed-forward transformation of expert e; $E$is the total number of experts; and HMoE(${\mathrm{h}}_i$) is the $Top-Kₑ$ weighted expert output.

Equation (11) renormalizes the gates over the selected experts, and Equation (12) computes the Top-2 weighted expert output. Experts outside ${\mathcal{T}}_i$ do not participate in the forward computation for token $i$. Sparse selection can, however, concentrate traffic on a small subset of experts. We therefore impose a balance constraint on the full routing probabilities before $Top-k$ truncation:

$${\overline{p}}_e^{\mathcal{l}}={\displaystyle \frac{1}{T}}{\sum }_{i=1}^T{p}_{i,e}^{\mathcal{l}}$$

(13)

$${\mathcal{L}}_{\mathrm{b}\mathrm{a}\mathrm{l}}^{\mathcal{l}}={\displaystyle \sum _{e=1}^E}{\left({\overline{p}}_e^{\mathcal{l}}−{\displaystyle \frac{1}{E}}\right)}^2$$

$${\mathcal{L}}_{\mathrm{b}\mathrm{a}\mathrm{l}}={\displaystyle \frac{1}{\left|{\mathcal{S}}_M\right|}}{\sum }_{\mathcal{l}\in {\mathcal{S}}_M}{\mathcal{L}}_{\mathrm{b}\mathrm{a}\mathrm{l}}^{\mathcal{l}}$$

(14)

where $T=BN$ is the number of node tokens in the current layer; i, e, and ℓ index tokens, experts, and HMoE layers, respectively; $p_{i,e}^{\mathcal{l}}$ is the soft routing probability before $Top-Kₑ$ truncation; ${\overline{p}}_e^{\mathcal{l}}$is the mean routing importance of expert e; $1/E$ is the uniform target importance; ${\mathcal{L}}_{\mathrm{b}\mathrm{a}\mathrm{l}}^{\mathcal{l}}$is the layer-wise balance loss; ${\mathcal{S}}_M$is the set of HMoE layers; $\left|{\mathcal{S}}_M\right|$ is its cardinality; and Lbal is the mean balance loss across those layers.

After HMoE aggregation, the L-layer backbone produces a final representation for each node. A lightweight $1\times 1$ regression convolution applies a shared linear projection to every node and uses $S$output channels to represent the $S$ future steps directly:

$${\widehat{Y}}_{b,:,n,1}=W_o{\mathrm{h}}_{b,n}^L+b_o\in {\mathbb{R}}^S,\widehat{Y}\in {\mathbb{R}}^{B\times S\times N\times 1}$$

(15)

where ${\mathrm{h}}_{b,n}^L$ is the final representation of node n in sample $b$; $W_o$ and $b_o$ are the shared regression weights and bias; the colon denotes all $S$ forecast horizons; the final index 1 denotes the single target channel; $\widehat{Y}$ is the prediction tensor; and $B,S$, and $N$ are the batch size, forecast length, and number of nodes, respectively.

During training, predictions are transformed back to the original traffic-flow scale. A masked mean absolute error is computed over valid targets and optimized jointly with the routing balance term:

$${\mathcal{L}}_{pred}={\displaystyle \frac{1}{\left|\mathsf{\Omega }\right|}}{\sum }_{(b,s,n)\in \mathsf{\Omega }}|{\widehat{Y}}_{b,s,n}-Y_{b,s,n}|,\mathcal{L}={\mathcal{L}}_{pred}+{\lambda }_{bal}{\mathcal{L}}_{bal},{\lambda }_{bal}=0.01$$

(16)

where $\mathsf{\Omega }$ is the set of valid target positions and $\left|\mathsf{\Omega }\right|$ is its cardinality; b, s, and n index the sample, forecast horizon, and sensor node; ${\widehat{Y}}_{b,s,n}$ and $Y_{b,s,n}$ are the prediction and ground truth; ${\mathcal{L}}_{pred}$ is the masked-MAE forecasting loss;${\lambda }_{bal}$ is the routing-balance loss; ${\lambda }_{bal}=0.01$ is its coefficient; and L is the total training objective. The first term directly optimizes forecast accuracy, while the second discourages persistent expert underuse. Direct multi-step output avoids recursive error accumulation and allows all horizons to share the same final node representation.

2.4. Datasets and Forecasting Protocol

We evaluate FAMoE-ST on the public PEMS04 and PEMS08 highway traffic-flow datasets derived from the California Performance Measurement System (PeMS). Both datasets are sampled at 5 min intervals and contain 307 and 170 sensor nodes, respectively. PEMS04 covers January–February 2018 and contains 16,992 time steps; PEMS08 covers July–August 2016 and contains 17,856 time steps. Table 1 summarizes the dataset statistics.

Samples are divided chronologically into training, validation, and test sets at a ratio of 6:2:2. For PEMS04, the three subsets contain 10,181, 3,393, and 3,395 samples; for PEMS08, they contain 10,699, 3,566, and 3,568 samples. Every experiment uses the preceding 12 steps to predict the next 12 steps; thus, 60 min of observations support direct forecasting over the subsequent 60 min. The traffic-flow channel is standardized using Z-score statistics computed exclusively from the training set. Predictions are transformed back to the original scale before evaluation, and positions with a ground-truth value of zero are excluded.

2.5. Evaluation Metrics and Implementation Details

Forecasting performance is evaluated using mean absolute error (MAE), root mean square error (RMSE), and mean absolute percentage error (MAPE), for which lower values indicate better performance. MAE summarizes absolute deviation, RMSE places greater weight on large errors, and MAPE measures relative error across traffic scales. Unless otherwise stated, the tables report the arithmetic mean of each metric over all 12 forecast steps.

The model uses a six-layer Transformer initialized from the first six blocks of a GPT-2 checkpoint. Because the sequence dimension represents sensor nodes, the native causal mask is replaced by an all-visible bidirectional attention mask. All node tokens use the shared zero GPT position ID, and learned node embeddings encode sensor identity without imposing an order-dependent spatial prior. The depth-dependent rule freezes attention and original FFN parameters in L1–L4, adapts L5–L6, and replaces only the FFNs in L5–L6 with HMoE layers. Each HMoE layer contains four experts and activates the Top-2 experts per token. The memory contains 32 learnable slots, of which the four most relevant are retrieved. The linear and memory routers are weighted equally (β = 0.5), and the balance-loss coefficient is 0.01. Input and output lengths are both 12, the batch size is 4, and models are trained for at most 200 epochs. We use the Ranger optimizer with an initial learning rate of 10⁻³, weight decay of 10⁻⁴, a gradient-clipping threshold of 5, and random seed 6666. The checkpoint with the lowest validation MAE is used for testing. For the capacity comparison, runtime is measured on an 8 GB NVIDIA GeForce RTX 4060 Ti with batch size 4 after five warm-up iterations and over 20 measured iterations. The environment uses Python 3.8.19, PyTorch 2.4.1+cu118, CUDA 11.8, cuDNN 9.1.0, and Transformers 4.30.2. All benchmarked operations use FP32 without automatic mixed precision. Each timed forward-backward or inference iteration is bracketed by torch.cuda.synchronize(), and peak allocated memory is obtained from torch.cuda.max_memory_allocated(). Mean epoch time is the wall-clock average over all 200 completed training epochs.

For consistent comparison, all baselines use the same chronological splits, 12-step input/output setting, masked evaluation metrics, and 200-epoch training budget, together with their recommended architectures and optimization settings. The baselines are Graph WaveNet (GWN) [26], Multivariate Time-Series Graph Neural Network (MTGNN) [27], Adaptive Graph Convolutional Recurrent Network (AGCRN) [28], Diffusion Convolutional Recurrent Neural Network (DCRNN) [29], Spectral Temporal Graph Neural Network (StemGNN) [30], Spatio-Temporal Graph Convolutional Network (STGCN) [31], Graph Multi-Attention Network (GMAN) [32], and Spatio-Temporal Normalization (STNorm) [33]. GATGPT and GCNGPT are graph-attention and graph-convolution variants of the ST-LLM framework and therefore follow the source implementation described in [10]. Table 2 reports MAE, RMSE, and MAPE at horizons 3, 6, and 12 and their 12-step averages.

2.6. Use of Generative Artificial Intelligence

During manuscript preparation, OpenAI Codex was used solely for Chinese-to-English translation and language editing. The authors reviewed and revised all generated text and assume full responsibility for the content of the manuscript.

3. Results and Discussion

3.1. Comparative Results

Table 2 compares the forecasting models at horizons 3, 6, and 12 and over the 12-step average. The average metrics characterize overall forecasting quality, whereas the horizon-specific values reveal how errors evolve as the prediction range increases.

On PEMS04, FAMoE-ST achieves 12-step average MAE, RMSE, and MAPE values of 18.55, 30.35, and 12.91%, respectively. Its MAE and RMSE rank first and are 1.49% and 1.81% lower than those of the corresponding second-best models, GMAN and GWN. Its MAPE ranks second and is only 0.22 percentage points above STNorm. FAMoE-ST therefore provides its most consistent advantage on PEMS04 in absolute-error control rather than in every relative-error comparison.

On PEMS08, FAMoE-ST obtains average MAE, RMSE, and MAPE values of 14.57, 24.08, and 9.71%, respectively. MAE and MAPE rank first: MAE is 1.62% lower than that of GMAN, while MAPE is 0.09 percentage points below the second-best GWN result. RMSE ranks second and is 0.39 higher than that of GWN. These results indicate improved average absolute accuracy while retaining competitive large-error sensitivity.

Across both datasets, FAMoE-ST outperforms the graph-enhanced GPT baseline GATGPT on all three average metrics and remains competitive at individual horizons. The clearest gains occur in MAE and RMSE, whereas MAPE does not dominate every baseline. Together with the controlled ablations in Section 3.3, these comparisons support the effectiveness of restricted upper-layer adaptation and hybrid-memory expert computation within the GPT-style backbone. They should not, however, be interpreted as evidence that linguistic pretraining itself causes the improvement.

To examine how absolute error evolves with forecasting range, Figure 2 plots horizon-wise MAE for FAMoE-ST and three representative baselines: AGCRN, DCRNN, and GCNGPT. These methods provide adaptive-graph recurrent, diffusion-recurrent, and graph-enhanced GPT references, respectively. This subset is used to maintain visual clarity rather than to present an exhaustive ranking; complete comparisons with all evaluated methods are reported in Table 2.

FAMoE-ST maintains the lowest MAE at every forecast horizon on both datasets. On PEMS04, its MAE increases from 17.05 at horizon 1 to 19.70 at horizon 12, whereas AGCRN, DCRNN, and GCNGPT reach 21.08, 22.57, and 24.26, respectively. The gap to the strongest displayed baseline therefore increases from 0.21 at horizon 1 to 1.37 at horizon 12. On PEMS08, FAMoE-ST increases from 12.81 to 15.76, compared with horizon-12 MAE values of 17.36, 18.11, and 20.39 for AGCRN, DCRNN, and GCNGPT. The corresponding gap to the strongest displayed baseline widens from 0.57 to 1.60. This consistent separation indicates that the proposed model accumulates error more slowly as the forecasting range grows.

Figure 3 provides a qualitative comparison of horizon-12 predictions over a continuous 24 h test segment. For each dataset, the displayed sensor is the node closest to the 75th percentile of ground-truth variability, and the non-overlapping 288-step segment has the highest ground-truth variability for that node. The local-view interval is likewise determined only from observed traffic dynamics; model predictions are not used in selecting the displayed case.

Across both examples, FAMoE-ST follows the principal rising and falling trends while reducing the lag or oversmoothing visible in the recurrent graph baselines. Under the same zero-masked evaluation protocol, local MAE on PEMS04 is 24.05 for FAMoE-ST, 25.31 for AGCRN, and 30.31 for DCRNN; the corresponding PEMS08 values are 16.32, 19.08, and 20.86. The enlarged intervals further show that FAMoE-ST responds more closely to rapid transitions and peak changes. These selected trajectories provide qualitative evidence and do not replace the full-test-set comparisons in Table 2.

3.2. Prediction Consistency and Data-Ratio Sensitivity

We additionally report the coefficient of determination (R²) to assess agreement between observed and predicted values; values closer to 1 indicate stronger agreement.

Figure 4a,b present the scatter relationships between observed and predicted flow in the training and test sets. Most points from both subsets lie near the y = x reference line, indicating close agreement across the principal range of traffic flow. On PEMS04, the training/test R² values are 0.966/0.960 and the corresponding RMSE values are 28.98/31.83; on PEMS08, they are 0.977/0.973 and 21.98/24.26. Figure 4 pools samples from all 12 horizons before computing RMSE, whereas Table 2 reports the arithmetic mean of 12 horizon-wise RMSE values; this difference in aggregation explains the small numerical discrepancy. From training to testing, R² decreases by only 0.006 and 0.004, while RMSE increases by 2.85 and 2.28. The observed train-test gaps are limited, although dispersion increases at high flows, identifying peak traffic as the principal remaining source of error.

Figure 4c,d show how R² and RMSE vary with the fraction of training data. On PEMS04, test R² ranges from approximately 0.959 to 0.963 and RMSE from 30.7 to 31.8. On PEMS08, R² decreases from approximately 0.978 to 0.973, while RMSE increases overall from about 22 to 24.3. Because performance is not monotonic in the data fraction, subset composition and stochastic training variation are plausible contributors. Nevertheless, R² remains above 0.95 at every ratio, indicating that the fitted relationship is retained under restricted training samples. Each curve is based on a single split; therefore, apparent local advantages at smaller ratios should be confirmed using multiple random seeds.

3.3. Ablation Study

To isolate the contributions of hierarchical freezing, sparse experts, memory routing, and balance loss, we compare the complete FAMoE-ST model with five variants: ST-LLM, ST-LLM without freezing, ST-LLM with linear MoE, FAMoE-ST without balance loss, and FAMoE-ST without freezing. Only the component combinations specified in Table 3 are changed; all remaining training and evaluation settings are held constant. Table 3 reports 12-step averages, and Figure 5 expresses each error change relative to the complete model.

The complete FAMoE-ST obtains MAE, RMSE, and MAPE values of 18.545, 30.353, and 12.914% on PEMS04 and 14.570, 24.083, and 9.710% on PEMS08. Removing frozen attention while retaining HMoE increases the three errors by 2.00%, 2.13%, and 0.26% on PEMS04 and by 2.52%, 1.25%, and 4.98% on PEMS08. Removing freezing from the single-FFN ST-LLM also worsens all three metrics on both datasets. Under the reported training protocol, restricting adaptation to the upper two layers therefore provides a consistent optimization benefit.

Adding linear sparse experts to ST-LLM reduces PEMS04 MAE from 19.935 to 19.534 and MAPE from 14.427% to 13.943%. On PEMS08, MAE decreases from 15.733 to 15.640 and MAPE from 10.414% to 10.218%, whereas RMSE changes only from 24.849 to 24.843. Adding memory routing without balance loss subsequently reduces MAE, RMSE, and MAPE by 4.67%, 2.37%, and 8.17% on PEMS04 and by 6.13%, 2.67%, and 4.44% on PEMS08 relative to linear MoE. Thus, linear MoE provides an initial gain over the shared FFN, while the memory branch accounts for most of the subsequent reduction under this ablation sequence.

The effect of balance loss is metric dependent. On PEMS08, adding the term reduces MAE, RMSE, and MAPE from 14.682, 24.180, and 9.764% to 14.570, 24.083, and 9.710%, corresponding to improvements of 0.76%, 0.40%, and 0.55%. On PEMS04, MAE and RMSE improve by 0.41% and 0.29%, but MAPE increases from 12.804% to 12.914%, or by 0.11 percentage points, consistent with the negative bar in Figure 5. Balance regularization therefore improves load allocation without guaranteeing a uniform gain for every dataset-metric pair.

Overall, the complete model reduces MAE, RMSE, and MAPE relative to ST-LLM by 6.97%, 3.18%, and 10.49% on PEMS04 and by 7.39%, 3.08%, and 6.76% on PEMS08. The sequence of ablations supports complementary contributions from restricted upper-layer adaptation, sparse expert capacity, and memory-assisted routing. Balance loss primarily improves absolute errors and the overall PEMS08 result, but its effect is not uniform across datasets and metrics.

3.4. Hybrid-Memory Routing Behavior

Figure 6 examines hybrid-memory routing from four complementary perspectives: expert usage across layers, token-routing partitions, the contributions of the linear and memory branches, and node-level expert preferences. Layers 1 and 2 in the figure correspond to backbone layers L5 and L6, respectively. Panels b and d analyze L5, whereas panel c analyzes L6. Because expert indices are permutation invariant, assignments are compared only within each dataset.

For each token in an input window, the model produces a fused routing distribution and a Top-1 expert assignment. The horizontal axes of the routing maps are reordered by dominant expert and preference strength solely to expose contiguous routing partitions; this ordering is not used at inference time. Branch contributions are averaged over tokens, whereas node preferences are aggregated over samples and time. Accordingly, Figure 6a summarizes global load, Figure 6b,d show local stability and sample-dependent switching, and Figure 6c compares the two routing branches.

Figure 6b,d indicate that routing combines relatively stable node specialization with dynamic switching across samples. Contiguous regions emerge after reordering, showing that some nodes maintain persistent expert preferences; nevertheless, the expert selected for a given node can change with the input traffic state. PEMS04 exhibits finer-grained alternation among experts, whereas PEMS08 is dominated by Expert 3. The router therefore reflects both persistent road-segment heterogeneity and sample-level flow variation, although the visualization alone does not assign a unique traffic meaning to any expert.

Figure 6c demonstrates that the memory branch contributes information not present in the instantaneous route. In the second PEMS04 module, the linear, memory, and fused probabilities for Expert 2 are 4.01%, 31.16%, and 17.58%. On PEMS08, the corresponding values for Expert 1 are 0.0307%, 30.97%, and 15.50%. The value 0.000307 is displayed as 0.02 for visibility, but the original value is retained in the analysis. At the same time, memory routing does not erase preferences expressed by the linear branch: for PEMS08 Expert 3, the linear probability is 51.24% and the fused probability remains 30.63%. These patterns are consistent with convex fusion of complementary instantaneous and memory-based routing signals.

Taken together, the diagnostics show nonuniform but noncollapsed conditional specialization. Node preferences provide a relatively stable component, sample-level switching accommodates short-term state variation, and the memory branch supplies cross-batch routing references. The evidence confirms that the hybrid-memory mechanism changes expert allocation, but it does not yet justify a unique traffic-semantic label for each expert. Analyses stratified by peak period, congestion stage, and road attributes are needed to establish such interpretations.

3.5. Capacity-Matched Performance and Efficiency

Because HMoE stores multiple expert copies, an improvement over a conventional FFN could otherwise be attributed to increased parameter capacity. We therefore construct a dense control in which the FFNs of L5 and L6 are widened until their stored parameter count matches that of HMoE. The two variants use the same GPT-2-initialized backbone, bidirectional spatial attention, shared zero position ID, L1–L4 freezing policy, data split, optimizer, and 200-epoch training budget. The dense control replaces each L5–L6 FFN with one expanded GPT-2 MLP whose intermediate width is 12,691; this gives 19,506,835 stored parameters per modified layer, compared with 19,507,328 for HMoE. Both variants are retrained as a matched pair on PEMS04; Table 4 reports the resulting accuracy and theoretical active computation.

For the MAC calculation, the Transformer hidden size is 768, each expert has an intermediate width of 3072, four experts are stored, two experts are activated, and the memory contains 32 slots. The dense count is 2 × 768 × 12,691 = 19,493,376 MACs per token per modified layer. The HMoE count is 2 × 2 × 768 × 3072 + 768 × 4 + 768² + 32 × 768 + 32 × 4 = 10,054,784 MACs. These expressions include the two FFN matrix multiplications, expert/router projections, memory-query projection, slot similarity, and slot-to-expert mapping. Softmax, Top-K indexing, tensor movement, and expert-dispatch overhead are excluded from the theoretical count and are instead captured by the empirical measurements in Table 5.

The two models differ by fewer than 0.001 million parameters and therefore provide an effectively capacity-matched comparison. Relative to the dense control, HMoE reduces MAE, RMSE, and MAPE by 2.34%, 3.27%, and 2.81%, respectively. At the same time, its theoretical FFN-and-router computation decreases from 19.493 to 10.055 million MACs per token in each modified layer, a reduction of 48.42%. These results indicate that the forecasting gain is not explained solely by a larger stored parameter count; sparse routing uses comparable capacity more effectively and activates substantially less computation per token.

The empirical measurements are consistent with the theoretical comparison. HMoE reduces forward-and-backward time by 19.02%, training peak memory by 16.71%, mean epoch time by 26.75%, and batch inference latency by 17.95%; its batch throughput increases by 21.87%. However, single-sample latency increases from 8.03 to 13.16 ms. This result identifies a practical boundary: sparse routing is advantageous when batched parallelism amortizes expert-dispatch overhead, whereas a dense FFN remains faster for isolated single-sample requests.

4. Conclusions

This study presented FAMoE-ST, a traffic-forecasting model that structurally adapts a GPT-style Transformer through ordering-neutral bidirectional spatial attention, hierarchical freezing, and hybrid-memory experts. Sensor histories are represented as node tokens; all tokens share the same GPT positional ID and attend to one another without a causal mask, while learned node embeddings retain sensor identity. The model freezes attention and original feed-forward parameters in L1–L4, adapts L5–L6, and combines instantaneous linear routing with memory routing. Top-2 gating then applies token-dependent transformations to heterogeneous traffic patterns.

FAMoE-ST achieves the best 12-step average MAE and RMSE and the second-best MAPE on PEMS04; on PEMS08, it achieves the best average MAE and MAPE and the second-best RMSE. Ablations indicate complementary gains from hierarchical freezing, sparse experts, and memory routing. Balance loss improves MAE and RMSE on both datasets and MAPE on PEMS08, but slightly degrades MAPE on PEMS04. Routing diagnostics further show that several experts remain active in both upper HMoE layers, stable node preferences coexist with sample-level switching, and the memory branch supplements the instantaneous router. A parameter-matched dense control further shows that HMoE improves all three forecasting metrics while reducing theoretical active computation, batch latency, and peak memory. These findings support the proposed architectural adaptation and routing mechanism rather than a causal claim about language-derived knowledge.

The evaluation is limited to two highway datasets, a fixed 12-step setting, and the reported backbone initialization. It therefore does not independently establish the contribution of GPT-2 language pretraining relative to a randomly initialized isomorphic Transformer, nor does it characterize variability across random seeds. Future work will evaluate additional regions, forecast horizons, external-event conditions, and matched initialization controls; report statistical significance across repeated runs; examine sensitivity to the freezing ratio, memory weight, balance coefficient, and expert capacity; and relate routing behavior to peak periods, congestion stages, and road attributes.