Paper deep dive
MercerFlow: Flow Matching in a Kernel-Induced Latent Space for Probabilistic Forecasting
Ilya Kuleshov, Egor Serov, Alexey Zaytsev
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Recent work has shown that probabilistic flow matching for time series forecasting benefits from a data-matched prior. The resulting prior introduces local correlations, which a sequential architecture usually absorbs: a recurrent neural network (RNN), a structured state-space model (S4), or a Transformer. However, such a backbone costs GPU memory and time per epoch. A cheaper alternative is MLP-based latent-space flow matching: embed the time series via an invertible map to a single latent vector and learn the flow there, so a tabular MLP can treat the series as a set of features. The relationship between the prior and the choice of linear latent map is understudied in conditional flow matching (CFM) forecasting, yet we found it strongly affects performance. Fixed transforms such as Fourier or discrete cosine (DCT) are only well-conditioned for Ornstein--Uhlenbeck priors, while a principal-component (PCA) map fit to the data is a strong but training-set-dependent reference sensitive to train--test shift. Instead, we propose to use the Mercer eigenbasis of the prior kernel: it diagonalises the centred covariance exactly, decouples from training data, and adapts to non-stationary and periodic priors. On five GluonTS benchmarks (ETTh1, ETTh2, Weather, Electricity, Traffic) under a shared protocol with TSFlow, the resulting MLP matches or beats it on CRPS at about $4.7\times$ less training memory and $3.5\times$--$4.4\times$ less time per epoch.
Tags
Links
- Source: https://arxiv.org/abs/2610.06039v1
- Canonical: https://arxiv.org/abs/2610.06039v1
Trouble viewing inline? Open PDF directly →
Full Text
105,152 characters extracted from source content.
Expand or collapse full text
MercerFlow: Flow Matching in a Kernel-Induced Latent Space for Probabilistic Forecasting Ilya Kuleshov †thanks: Equal contribution.†thanks: i.kuleshov@applied-ai.ru Egor Serov11footnotemark: 1 Alexey Zaytsev Affiliation: Applied AI Institute Affiliation: Moscow, Russia Abstract Recent work has shown that probabilistic flow matching for time series forecasting benefits from a data-matched prior. The resulting prior introduces local correlations, which a sequential architecture usually absorbs: a recurrent neural network (RNN), a structured state-space model (S4), or a Transformer. However, such a backbone costs GPU memory and time per epoch. A cheaper alternative is MLP-based latent-space flow matching: embed the time series via an invertible map to a single latent vector and learn the flow there, so a tabular MLP can treat the series as a set of features. The relationship between the prior and the choice of linear latent map is understudied in conditional flow matching (CFM) forecasting, yet we found it strongly affects performance. Fixed transforms such as Fourier or discrete cosine (DCT) are only well-conditioned for Ornstein–Uhlenbeck priors, while a principal-component (PCA) map fit to the data is a strong but training-set-dependent reference sensitive to train–test shift. Instead, we propose to use the Mercer eigenbasis of the prior kernel: it diagonalises the centred covariance exactly, decouples from training data, and adapts to non-stationary and periodic priors. On five GluonTS benchmarks (ETTh1, ETTh2, Weather, Electricity, Traffic) under a shared protocol with TSFlow, the resulting MLP matches or beats it on CRPS at about 4.7×4.7× less training memory and 3.5×3.5×–4.4×4.4× less time per epoch. 1 Introduction The usual approach to probabilistic forecasting is to devise an algorithm which samples from a predictive distribution over future windows given the observed past. A sampler that is expressive, while staying cheap to train and cheap to run, is the practical target. Flow matching (Lipman et al., 2023; Tong et al., 2024) is a prominent modern way to build such samplers: one learns a velocity field along a probability path from a tractable prior to the data, and training is ordinary regression onto a known interpolant. Since each forecast requires numerically integrating this velocity during inference, the cost at test time is dominated by the backbone that evaluates the field; cheapness at training time matters too, since the backbone also dictates GPU memory and time per epoch. If the prior already matches the data, the flow starts closer to the target and needs less capacity to finish the trip. TSFlow (Kollovieh et al., 2025) commits to a Gaussian-process (GP) prior on the forecast window for that reason. However, under that prior, neighbouring time coordinates stay correlated, so a global multilayer perceptron (MLP) spends capacity on local smoothness the GP already encodes. TSFlow absorbs those correlations with a bidirectional state-space backbone (Gu et al., 2022), which is accurate on long smooth series but expensive in memory and time per epoch. Orthonormal projections are a classical way to decorrelate coordinates in a latent space (Mercer, 1909; Karhunen, 1947; Gabor, 1946). Several recent generative models run denoising or score matching in orthonormal latents: DCTdiff (Ning et al., 2025) trains image diffusion in discrete cosine coordinates, Crabbé et al. (Crabbé et al., 2024) dualise time-series diffusion to the Fourier domain, and Phillips et al. (Phillips et al., 2022) score functional data along truncated Mercer eigenfunctions. However, when we commit to a GP prior, as in TSFlow, the latent map must also decorrelate the resulting structure—not just the correlations in the data itself. To address this, we conduct a controlled comparison of five linear codecs for projection on a shared architecture and matched prior: Fourier (discrete Fourier transform), discrete cosine (DCT), principal components (PCA), the identity (Raw), and Mercer-based eigenbasis. Mercer’s theorem supplies the orthonormal eigenbasis of the kernel one has already committed to (Mercer, 1909; Karhunen, 1947), so the coefficients are uncorrelated under the corresponding GP. Thus, that decorrelation is exact for the centred prior and approximate for the conditional law and for the interpolant on which the network trains. Across codecs, Mercer is the most versatile data-free choice: it works on data which requires a non-OU prior (Fourier and DCT best diagonalise OU-like priors, which are not always optimal for the data) and where PCA suffers from the train–test shift since it fits only the training path (§3.2). That said, all five codecs are orthonormal projections, so a rotation-equivariant optimiser with an MSE loss should treat them as equivalent reparameterisations (Nesterov, 2004; Zhang et al., 2025). However, that equivalence is optimiser-specific: it holds for Euclidean SGD (Proposition 1), but fails under Adam because of the first-layer input geometry. Write A=[⊤]A=E[xx ] for a first-layer input x second-moment matrix and κ(A)=λmax(A)/λmin(A)κ(A)= _ (A)/ _ (A) for the condition number of its spectrum. Orthonormal codecs leave the spectrum and κ(A)κ(A) unchanged, so Euclidean SGD cannot prefer one codec over another. Under Adam, however, per-coordinate rescaling is modelled by the Jacobi preconditioner D=diag(A)−1/2D=diag(A)^-1/2, after which the normalised second-moment matrix is DADDAD with the geometry conditioned by κ(DAD)κ(DAD). The rotation leaves κ(A)κ(A) fixed but changes diag(A)diag(A), which is what Jacobi reads, so κ(DAD)κ(DAD) is codec-dependent. Empirically, the decorrelating codecs exhibit both lower κ(DAD)κ(DAD) and better forecasting accuracy than Raw under Adam on challenging datasets, while their advantage disappears under SGD. Our experiments therefore examine how codec choice interacts with the optimiser and whether MercerFlow offers a favourable balance of predictive accuracy and computational cost relative to TSFlow and other baselines. Figure 1: Overview of MercerFlow. The joint context and horizon window is mapped into the decorrelated Mercer eigenbasis of the GP prior. A lightweight MLP learns conditional flow matching directly in this latent space, enabling fast and scalable probabilistic forecasting without sequential backbones. Contributions. Our contributions are as follows: • A controlled five-codec study of conditional flow matching under matched architecture, prior, and protocol: Mercer, Fourier, DCT, PCA, and Raw on five GluonTS benchmarks with paired t-test ranks. • A diagnostic of codec gaps under Adam versus SGD. First-layer κ(DAD)κ(DAD) separates Raw from the decorrelating codecs on hard sets under Adam, while the gaps vanish when the models are trained with SGD. • MercerFlow: conditional flow matching carried out directly in the Mercer / Karhunen–Loève eigenbasis of the forecasting GP prior. Unlike fixed transforms such as DCT and Fourier, the basis adapts directly to non-stationary and periodic priors, providing a data-free representation that remains robust across kernels without fitting the training path. • An empirical comparison against TSFlow on multiple benchmarks. The MLP codecs are about 4.7×4.7× cheaper in peak GPU memory and 3.5×3.5×–4.4×4.4× cheaper in time per epoch at matched batch size, making MLP flow matching in linearly projected latents a practical alternative to sequential S4 backbones. 2 Method We study probabilistic time series forecasting, aiming to model the conditional distribution of a future window given the observed past. Following the channel-independent setting of TSFlow (Kollovieh et al., 2025), each channel is treated as a separate univariate series. Let LpL_p and LfL_f denote the context (past) and forecast horizon (future) lengths, respectively, yielding a total sequence length L=Lp+LfL=L_p+L_f. We define the regular discrete time grid =0,…,L−1S=\0,…,L-1\, partitioned into the observed past p=0,…,Lp−1S_p=\0,…,L_p-1\ and the future forecast horizon f=Lp,…,L−1S_f=\L_p,…,L-1\. The observed past series is denoted by p∈ℝLpy_p ^L_p, the true future window by f∈ℝLfy_f ^L_f, and their joint concatenation by =[p,f]∈ℝLy=[y_p,y_f] ^L. MercerFlow overview. Figure 1 summarises three stages of MercerFlow. First, we condition a Gaussian-process prior on the observed past and sample an initial future window (Subsection 2.1). Second, we concatenate this sample with the past and encode the joint window in the orthonormal eigenbasis of the same kernel matrix (Subsection 2.2). Third, during training, a lightweight MLP learns the velocity along linear interpolants between the encoded source and target windows (Subsection 2.3). Finally, at inference time, we integrate the learned velocity from a newly sampled source, apply the inverse transform, and extract the forecast horizon (Subsection 2.3). Repeating this procedure produces samples from the predictive distribution. 2.1 Conditional GP Prior Let k:×→ℝk:S×S be a positive-definite covariance kernel, e.g., Ornstein–Uhlenbeck (OU), radial basis function (RBF), or periodic, with length scale ℓ>0 >0 and diagonal jitter ε>0 >0. In our primary benchmark runs, we adopt the stationary exponential (OU) kernel: k(s,s′)=exp(−|s−s′|ℓ),s,s′∈.k(s,s )= (- |s-s | ), s,s . (1) Evaluating k across all pairs of temporal coordinates in S defines the joint covariance matrix K∈ℝL×LK ^L× L, with entries Kij=k(i,j)+εI[i=j]K_ij=k(i,j)+ I[i=j]. Partitioning K according to the past and future index sets pS_p and fS_f yields: K=[KppKpfKfpKff],K= bmatrixK_p&K_pf\\ K_fp&K_f bmatrix, (2) where Kpp∈ℝLp×LpK_p ^L_p× L_p, Kff∈ℝLf×LfK_f ^L_f× L_f, and Kfp=Kpf⊤∈ℝLf×LpK_fp=K_pf ^L_f× L_p. Conditioned on the observed context py_p, the observation-space prior over the future window 0∈ℝLfy_0 ^L_f is given by the standard Gaussian Process posterior: 0∼(,),=KfpKpp−1p,=Kff−KfpKpp−1Kpf.y_0 ( μ, ), μ=K_fpK_p^-1y_p, =K_f-K_fpK_p^-1K_pf. (3) We select the prior length scale ℓ via validation-likelihood maximisation per dataset (Table 8, §3). The diagonal jitter ε provides numerical stability but cannot substitute for basis rotation: as shown in Appendix F.3, Raw’s condition number remains near 10510^5 across all ε≤10−2 ≤ 10^-2, while extreme regularisation (ε=1 =1) degenerates the prior into white noise and degrades both Mercer and Raw equally. 2.2 Linear Codec on the Joint Window To make the observed past available to the velocity field, we apply a single linear codec to the joint context and future window. The source window concatenates py_p with a sample 0y_0 from the conditional GP prior, while the target window combines the same context with the observed future fy_f. Both windows are projected into a shared coefficient space using an orthonormal matrix Φ∈ℝL×L ^L× L, with encoding and decoding defined by: =Φ⊤,=Φ.c= y, = . (4) Since the columns of Φ are orthonormal, ΦΦ⊤=IL =I_L, and this map is an exact isometry. The source prior state is 0=Φ⊤[p,0]c_0= [y_p,y_0] and the target state is 1=Φ⊤[p,f]c_1= [y_p,y_f]. The five codecs described below differ only in how they construct the basis vectors ϕi φ_i of the matrix Φ=[ϕ1,…,ϕL] =[ φ_1,…, φ_L]. Identity (Raw). This codec keeps the series in observation coordinates with Φ=IL =I_L, which yields =c=y. Each coordinate represents an ambient time step. It acts as an unrotated baseline and matches the time-domain setup of TSFlow, but uses an MLP rather than an S4 backbone. Mercer (Ours). The basis columns ϕi φ_i are the orthonormal eigenvectors of the kernel matrix K: K=ΦΛΦ⊤,Λ=diag(λ1,…,λL).K= , =diag( _1,…, _L). (5) Eigenvalues are sorted in descending order. The first row of Φ is constrained to be non-negative to fix signs. By Mercer’s theorem and the Karhunen–Loève expansion, if ∼(,K)y (0,K), then Cov()=ΛCov(c)= (Rasmussen and Williams, 2006). Appendix B.3 gives the exact derivation. Therefore, latent coordinates are uncorrelated under the prior. This decorrelation is exact for the prior and provides an effective preconditioning proxy for the conditional law. Fourier. This codec implements a real discrete Fourier transform on the time grid s∈s (Gabor, 1946; Baraniuk and Jones, 1995). The first column is the constant direct current (DC) mode L−1/2L^-1/2. Higher harmonic frequencies m=1,2,…m=1,2,… provide cosine and sine pairs: ϕmcos(s)=2Lcos(2πmsL),ϕmsin(s)=2Lsin(2πmsL). _m (s)= 2L\, ( 2π msL ), _m (s)= 2L\, ( 2π msL ). (6) When L is even, the final column is the Nyquist mode L−1/2cos(πs)L^-1/2 (π s). The transform =Φ⊤c= y expresses the series in this real orthonormal Fourier basis, with each coefficient ci=ϕi⊤c_i= φ_i y corresponding to one Fourier mode. DCT. This codec uses an orthonormal type-I discrete cosine transform, ordered from lowest to highest frequency. Evaluated at time steps s∈s , the basis columns are: ϕ0DCT(s)=L−1/2,ϕmDCT(s)=2Lcos(πm(s+1/2)L),m=1,…,L−1. _0^DCT(s)=L^-1/2, _m^DCT(s)= 2L\, ( π m(s+1/2)L ), m=1,…,L-1. (7) With L−1/2L^-1/2 and 2/L 2/L normalisation, the modes form an orthonormal basis. Unlike Fourier, the DCT does not enforce periodic boundaries. It provides an asymptotic Karhunen–Loève basis for stationary Markovian and OU processes (Ahmed et al., 1974; Gray, 2006). PCA. The basis columns ϕi φ_i are the eigenvectors of the empirical sample covariance matrix computed across all training windows (k)k=1S\y^(k)\_k=1^S (Jolliffe, 2005): Σ^train=1S∑k=1S((k)−¯)((k)−¯)⊤,Σ^train=ΦΛΦ⊤, _train= 1S _k=1^S(y^(k)- y)(y^(k)- y) , _train= , (8) here ¯ y is the empirical mean of the training sequences. Eigenvalues are sorted in descending order with the same sign convention as Mercer. Centering is used solely to construct Σ^train _train, while the latent projection =Φ⊤c= y remains strictly linear. Because Φ is fitted to the training sample, PCA is data-dependent and sensitive to train–test distribution shifts. 2.3 Flow Matching in Coefficient Space Flow matching (Lipman et al., 2023) learns a time-dependent velocity field that transports a source distribution to the data. Standard formulations often require minibatch optimal transport to straighten trajectories originating from unstructured Gaussian noise (Tong et al., 2024). In contrast, we adopt TSFlow (Kollovieh et al., 2025) and use the conditional GP prior equation 3 as the source distribution. Because this prior already captures temporal history and local smoothness, independent pairing between source and target yields straight paths without optimal transport overhead. We instantiate this probability path directly in the joint coefficient space: t=(1−t)0+t1,=1−0,c_t=(1-t)c_0+tc_1, =c_1-c_0, (9) with flow time t∼[0,1]t [0,1], source 0c_0, and target 1c_1. Since the historical context is shared, the target displacement simplifies to =Φ⊤[,f−0]u= [0,y_f-y_0]. We train an MLP velocity field θ(t,t)v_θ(t,c_t) by minimizing the mean squared regression error: ℒ(θ)=t,0,1‖θ(t,t)−‖22.L(θ)=E_t,c_0,c_1 \|v_θ(t,c_t)-u \|_2^2. (10) Because Φ is an orthogonal isometry (ΦΦ⊤=IL =I_L), it preserves the Euclidean norm: ‖2=‖f−0‖2\|u\|_2=\|y_f-y_0\|_2. Consequently, the training loss magnitude is identical across all five codecs (Appendix B.1). To preserve input-layer geometry, the MLP injects AdaLN time modulation only after the first linear map (Appendix C.3). During inference, we sample 0∼(,)y_0 ( μ, ), project the joint state 0c_0, and integrate θv_θ from t=0t=0 to 11 using an Euler ODE solver to obtain the terminal state ^1 c_1. The predicted horizon ^f y_f is then extracted from the reconstructed sequence ^=Φ^1 y= c_1. 2.4 Why the Rotation Can Matter: Adam’s Diagonal At full order, corresponding to a square orthonormal basis Φ∈ℝL×L ^L× L on the ambient sequence window, Mercer, Fourier, and the identity baseline parameterise the identical Euclidean interpolant in three distinct orthonormal coordinate systems. Proposition 1 (SGD is invariant to the choice of orthonormal codec). Let Φ∈ℝL×L ^L× L be a square orthogonal matrix defined on the ambient observation space. Consider the velocity field θv_θ parameterized by θ=W1,1,W2,2,θψθ=\W_1,b_1,W_2,b_2, _ψ\: θ(t,)=W2ψ(W1+1,t,θψ)+2,v_θ(t,z)=W_2\,ψ (W_1z+b_1,\,t;\, _ψ )+b_2, (11) where ψ(⋅,t,θψ)ψ(·,t;\, _ψ) aggregates all internal hidden layers along with AdaLN time conditioning, acting entrywise on its first argument with parameters θψ _ψ independent of Φ . Let θI=W1,1,W2,2,θψ _I=\W_1,b_1,W_2,b_2, _ψ\ denote the parameters of the identity network, and define the conjugated parameter set θΦ=W1Φ,1Φ,W2Φ,2Φ,θψΦ _ =\W_1 ,b_1 ,W_2 ,b_2 , _ψ \ of the codec-Φ network by W1Φ=W1Φ,W2Φ=Φ⊤W2,2Φ=Φ⊤2,1Φ=1,θψΦ=θψ.W_1 =W_1 , W_2 = W_2, _2 = b_2, _1 =b_1, _ψ = _ψ. (12) Then the velocity field satisfies θΦ(t,Φ⊤)=Φ⊤θI(t,)v_ _ (t, z)= v_ _I(t,z), and the squared conditional flow matching loss equation 10 is preserved identically: ‖θΦ−‖22=‖θI−(1−0)‖22,\|v_ _ -u\|_2^2=\|v_ _I-(z_1-z_0)\|_2^2, (13) where =Φ⊤(1−0)u= (z_1-z_0). Under rotation-invariant parameter initialisations, Euclidean gradient descent and Euclidean momentum schemes, including the heavy-ball method and Nesterov accelerated gradient, generate identical loss trajectories and identical decoded flows at every training iteration (Nesterov, 2004; Zhang et al., 2025). The proof is provided in Appendix B.1. An empirical validation using Nesterov SGD at a matched computational budget is reported for ETTh1 and ETTh2 in Table 3 of Appendix B.2. The equivalence established in Proposition 1 is strictly Euclidean, as standard stochastic gradient descent applies isotropic updates across all coordinate axes. Modern probabilistic forecasting models are instead routinely optimized with adaptive algorithms such as AdamW (Kingma and Ba, 2015). These optimizers maintain an exponential moving average of squared coordinate-wise gradients to rescale parameter updates individually, thereby dampening oscillations along steep directions and accelerating progress along plateaus. Because this coordinate-wise scaling is defined relative to the canonical parameter axes, it is not rotation-equivariant. Rotating the input-space basis fundamentally alters Adam’s trajectory, whereas SGD’s trajectory undergoes a simple linear conjugation (Zhang et al., 2025; Xie et al., 2025). Consequently, when training the velocity field θv_θ of Section 2.3 with AdamW, optimisation dynamics are governed not by the raw input Gram matrix, but by its geometry following diagonal preconditioning. Let ∈ℝLx ^L denote the input vector to the first linear layer, and let A=[⊤]A=E[xx ] be its uncentred second-moment matrix. For the linear mapping =Wh=Wx, the matrix A determines the curvature of the corresponding linear subproblem. The extremal directions of stretch are given by the eigenvectors of A, and their relative ratio is quantified by the condition number κ(A)=λmax(A)λmin(A).κ(A)= _ (A) _ (A). (14) As Nesterov (2004) show, gradient descent on a quadratic objective has a convergence rate governed directly by κ(A)κ(A). Because an orthogonal transformation leaves the spectrum of A invariant, substituting A with Q⊤AQ AQ preserves κ(A)κ(A) and cannot accelerate Euclidean gradient descent. In contrast, the first-layer dynamics of Adam deviate substantially from Euclidean optimisation. Under the Kronecker-factored curvature model of Martens and Grosse (2015), the coordinate-wise second-moment accumulator decomposes into an output factor and an input factor given by diag(A)diag(A). The resulting coordinate rescaling corresponds to the classical Jacobi preconditioner D=diag(A)−1/2D=diag(A)^-1/2, producing the preconditioned input ~=D x=Dx and the rescaled curvature matrix DADDAD. The optimisation of the first-layer weights is therefore dictated by κ(DAD)κ(DAD). According to van der Sluis (1969), the Jacobi preconditioner minimises the condition number among all diagonal scalings up to a factor depending only on dimensionality. When A possesses small off-diagonal elements, DADDAD approximates the identity matrix, effectively whitening the layer inputs. Proposition 2 (Jacobi preconditioning distinguishes linear codecs). Let the total window length be n=Ln=L, and let Tn(f)∈ℝn×nT_n(f) ^n× n denote the symmetric Toeplitz covariance matrix generated by the spectral density function f(ω)f(ω) of a stationary Ornstein–Uhlenbeck process with length scale ℓ : [Tn(f)]jk=12π∫−πf(ω)ei(j−k)ωω=exp(−|j−k|ℓ). [T_n(f) ]_jk= 12π _-π^πf(ω)\,e^i(j-k)ω\,dω= (- |j-k| ). (15) Let =Φ⊤x= y be the coefficient representation under the centred prior ∼(,K)y (0,K) with K=Tn(f)+δInK=T_n(f)+δ I_n, and let A=[⊤]A=E[xx ]. Under the Jacobi preconditioner D=diag(A)−1/2D=diag(A)^-1/2, the condition number κ(DAD)κ(DAD) satisfies: 1. Identity (Raw): The matrix A=Tn(f)+δInA=T_n(f)+δ I_n is Toeplitz with constant diagonal 1+δ1+δ. The preconditioner satisfies D=(1+δ)−1/2InD=(1+δ)^-1/2I_n, yielding κ(DAD)=κ(A)κ(DAD)=κ(A), which grows monotonically with the prior length scale ℓ . 2. Mercer: Because the columns of Φ diagonalize the kernel matrix K, the second-moment matrix satisfies A=ΛA= , where Λ is diagonal. Hence DAD=InDAD=I_n, which yields the optimal condition number κ(DAD)=1κ(DAD)=1. 3. Fourier and DCT: The respective transformation bases satisfy the asymptotic bounds κ(DAD)≤1+(1minωf(ω)),κ(DAD)≤ 1+O ( 1 _ωf(ω) ), (16) where the bound loosens whenever f(ω)f(ω) concentrates mass predominantly at low frequencies. The proof is given in Appendix B.3, where Table 4 provides a structured overview of the second-moment structure and condition numbers across all five codecs. On synthetic stationary OU paths where the data-generating kernel is known exactly, empirical κ(DAD)κ(DAD) matches this ordering and CRPS ties across codecs within 0.2%0.2\% (Appendix F.1), reproducing under analytical control the behavior observed on easy real datasets. While an orthonormal change of basis cannot alter κ(A)κ(A), it modifies the diagonal entries diag(A)diag(A) that govern Adam’s scaling factors. During training, the network input is not drawn from the centred prior y, but represents the conditional interpolant tc_t. Empirical evaluations of κ(DAD)κ(DAD) across datasets are detailed in Appendix C.5 (Table 12) and summarized in Section 3.2. Under AdamW optimisation, first-layer κ(DAD)κ(DAD) clearly separates the ill-conditioned Raw baseline from the decorrelating codecs on challenging benchmarks, whereas training with SGD eliminates performance discrepancies across codecs, in complete alignment with Proposition 1. 3 Experiments Benchmark Setup and Protocol. Experiments cover five standard time series benchmarks: ETTh1, ETTh2, Weather, Electricity, and Traffic. The context and forecast lengths are fixed at Lp=Lf=192L_p=L_f=192, giving a joint window of length L=384L=384, with channels treated independently. Baselines include four alternative codecs (Raw, Fourier, DCT, and PCA), the flow-based forecaster TSFlow, and the deterministic forecasters DLinear, PatchTST, and DMamba. All five codec variants use the same MLP architecture, joint-window conditioning scheme, and GP prior. Prior hyperparameters are selected by validation likelihood (Appendix C.4). Architecture and training details appear in Appendix C. For each probabilistic model, evaluation uses S=100S=100 forecast trajectories ^f(s)s=1S\ y_f^(s)\_s=1^S per test window. Normalised continuous ranked probability score (CRPS) measures distributional forecast accuracy (Alexandrov et al., 2020; Kollovieh et al., 2025), while mean absolute error (MAE) measures point-forecast accuracy. For probabilistic models, the point forecast is the empirical median of the generated samples. Scores are aggregated over rolling test windows, with metric definitions and aggregation details provided in Appendix C.1. 3.1 Main Results: Forecasting Quality and Efficiency Table 1 summarises overall prediction quality across all five benchmarks via paired two-sided t-test ranks, while Table 2 reports the computational footprint. Full per-dataset numerical CRPS and MAE distributions are reported in Appendix D within Tables 14 and 15. Table 1: Predictive performance summary across all five benchmarks (ETTh1, ETTh2, Weather, Electricity, Traffic). Rank =1+=1+ number of methods strictly outperforming the model under a paired two-sided t-test on shared seeds (p<0.2p<0.2, lower score is better); a win denotes rank 11. Ours Linear Codec Baselines Sequential Metric Mercer Fourier DCT PCA Raw TSFlow (S4) Mean CRPS rank (↓ ) 1.001.00 1.201.20 1.401.40 1.801.80 3.403.40 5.805.80 # CRPS wins (↑ ) /5/5 4/54/5 4/54/5 3/53/5 2/52/5 0/50/5 Mean MAE rank (↓ ) 1.001.00 1.601.60 1.401.40 2.002.00 3.203.20 5.005.00 # MAE wins (↑ ) /5/5 3/53/5 4/54/5 3/53/5 2/52/5 1/51/5 Table 2: Computational efficiency benchmark at batch size 512. Comparison of per-epoch training time (seconds, ↓ ) and peak GPU memory (GB, ↓ ). Metric / Architecture ETTh1 ETTh2 Weather Electricity Traffic MLP Time (s, ↓ ) 0.67 0.69 1.80 0.93 0.68 S4 Time (s, ↓ ) 2.50 2.40 7.70 4.10 2.60 MLP VRAM (GB, ↓ ) 1.5 1.5 1.5 1.5 1.5 S4 VRAM (GB, ↓ ) 7.1 7.1 7.1 7.1 7.1 Probabilistic Forecasting Quality. As shown in Appendix D, flow matching models decisively outperform deterministic architectures on distributional calibration. While deterministic baselines achieve competitive point estimates on simple trends, they degrade on multi-frequency series (e.g., on Traffic, MercerFlow achieves 8.70×10−28.70× 10^-2 CRPS versus 10.7510.75 for DLinear and 10.4210.42 for DMamba; on Weather, 1.821.82 versus 6.266.26). Compared to TSFlow, which employs a heavy bidirectional S4 backbone, our lightweight MLP codecs achieve substantially better CRPS calibration across all benchmarks (mean rank 1.001.00–1.801.80 vs 5.805.80, Table 1). By offloading temporal continuity directly to the analytical GP prior, an MLP operating in decorrelated coordinates avoids wasting capacity on local correlations. Within the MLP family, the unrotated Raw baseline consistently degrades on harder datasets (ETTh1, Electricity, Traffic), yielding an inferior mean CRPS rank of 3.403.40. In contrast, MercerFlow resolves first-layer conditioning bottlenecks, securing a flawless mean rank of 1.001.00 (5/55/5 wins) across both CRPS and MAE. The full per-dataset CRPS and MAE scores, including the deterministic baselines and the horizon extension to Lf∈336,720L_f∈\336,720\, appear in Appendix D (Tables 14, 15, and 16). Computational Efficiency. As detailed in Table 2, operating in the latent coefficient space allows a compact feed-forward MLP (≈55k≈ 55k parameters) to replace the sequential state-space model (186k186k parameters). This achieves a 4.7×4.7× reduction in peak GPU memory (1.51.5 GB vs 7.17.1 GB) and accelerates per-epoch throughput by 3.5×3.5×–4.4×4.4× at matched batch size (512512). Furthermore, straight probability paths allow MercerFlow to saturate predictive performance within only 4–10 ODE integration steps, drastically outperforming equivalent score-based diffusion models (Appendix F.2, Table 19). 3.2 Why Mercer? Codec Universality and Prior Alignment Traffic uniquely exhibits strong weekly periodicity: validation GP likelihood decisively selects a weekly periodic kernel (p=168p=168, ℓ=0.1 =0.1, NLL 0.0350.035) over the stationary OU prior (ℓ=7 =7, NLL 0.4650.465; Table 6). We therefore use Traffic to compare the default OU prior against a better periodic prior, exposing the differences between fixed harmonic transforms, empirical bases, and MercerFlow’s analytical basis. Figure 2 shows first-layer Jacobi Gram heatmaps on conditional-prior draws for both Traffic priors. Under the stationary OU prior, DCT is near-optimal because its fixed harmonic basis approximates the Karhunen–Loève expansion of Markovian processes. Under the periodic kernel, DCT cannot adapt; off-diagonal energy concentrates at weekly harmonics, surging DCT’s condition number κ(DAD)κ(DAD) by 36×36× (1.4×1031.4× 10^3 for Mercer, 5.2×1045.2× 10^4 for DCT) and dropping its test CRPS (Table 10). MercerFlow recomputes Φ analytically from the chosen kernel K, diagonalising the periodic kernel by construction. While Fourier avoids DCT’s mirror reflection mismatch on periodic series, it suffers from severe boundary-leakage ill-conditioning on stationary long-horizon series (κ(DAD)κ(DAD) reaching 268.5268.5 on ETTh1 vs 1.001.00 for Mercer) and becomes unstable under learning rate stress on the periodic prior (Table 10). While PCA attains low training condition numbers, it overfits empirical samples and degrades on held-out test windows (Table 12: PCA condition number surges from 481481 to 12901290, whereas Mercer shifts only from 683683 to 925925). Appendix C.5 details full κ(DAD)κ(DAD) tables across all benchmarks. Figure 2: Traffic decorrelation on conditional-prior draws 0=Φ⊤[p;0]c_0= [y_p;y_0] for the joint window. Each panel displays |DAD||DAD| for the leading 6464 modes, with condition number κ(DAD)κ(DAD) and test CRPS (η=10−3η=10^-3) shown above each plot. Aligning the codec basis with the prior kernel maintains low off-diagonal correlation and preserves downstream performance. 3.3 Extended Analyses and Appendix Overview We conduct extensive controlled experiments and ablations in the appendix: Optimizer Equivariance under SGD. Proposition 1 predicts Euclidean SGD is rotation-equivariant; matched-budget Nesterov SGD ties all five codecs within one standard deviation on ETTh1 and ETTh2 (Appendix B.2, Table 3). Under AdamW, coordinate-wise rescaling breaks rotational equivariance, unlocking large gains for well-conditioned codecs on hard benchmarks. High Learning Rate Dynamics and Stability. Aggressive learning rates (η=10−1η=10^-1) collapse the unrotated baseline (CRPS explodes to 1.221.22 on ETTh1 and 0.2460.246 on Traffic; Table 17), whereas all decorrelating codecs remain stable (Appendix E, Figure 3). Synthetic OU and Diffusion Comparison. On synthetic OU paths where the true kernel is known, Mercer achieves κ(DAD)=1.00κ(DAD)=1.00, and misspecifying the length scale by 10×10× degrades CRPS by only 0.0060.006 (Appendix F.1, Table 18). Furthermore, straight probability paths saturate calibration in 44–1010 integration steps, outperforming score-based diffusion at 128128 steps (Appendix F.2, Table 19). Prior Regularisation and Full Scoreboard. Diagonal jitter δ∈[10−4,1.0]δ∈[10^-4,1.0] confirms that damping the prior cannot rescue an unrotated basis: Raw stays ill-conditioned (κ≈105κ≈ 10^5; Appendix F.3, Table 20). Full per-dataset CRPS/MAE scoreboards and extended horizons (Lf∈336,720L_f∈\336,720\) appear in Appendix D (Tables 14, 15, and 16). 4 Conclusion Gaussian process priors provide structured probability paths for conditional flow matching, but raw temporal coordinates leave dense correlations that create severe conditioning bottlenecks under AdamW. First-layer preconditioning separates these representations on challenging benchmarks: the unrotated baseline consistently trails behind and collapses under learning rate stress, while rotation-equivariant SGD leaves all orthonormal codecs tied. Fixed harmonic transforms fail when the covariance departs from stationary Toeplitz structures, and empirical PCA degrades under distribution shifts between train and test windows. By contrast, the Mercer eigenbasis diagonalizes the chosen prior covariance analytically without data fitting, maintaining stable preconditioning across both stationary and periodic processes. Aligning the latent coordinate frame with the prior kernel allows a lightweight feed-forward MLP to match or surpass sequential state-space models on probabilistic calibration, while reducing peak GPU memory by 4.7×4.7× and per-epoch training time by 3.5×3.5×–4.4×4.4×. Limitations Our theoretical analysis idealises Adam’s coordinate rescaling as exact Jacobi preconditioning and assumes rotation-invariant SGD initialisation. Empirically, performance depends on selecting an appropriate prior kernel and length scale via validation likelihood; also, our evaluations are restricted to univariate channels on a single fixed-width MLP backbone. Finally, the condition number metric κ(DAD)κ(DAD) captures the relative performance among Mercer, DCT, Fourier, and PCA bases primarily under structured departures from stationary Toeplitz priors (such as elevated learning rates, periodic kernels, or train–test distribution shifts). This aligns with our theoretical analysis, which shows that preconditioning differences diminish when the kernel matches the data-generating process. Reproducibility We are publishing our entire repository on an anonymous repository to facilitate reproducibility: https://anonymous.4open.science/r/mercerflow_public-0323/. All datasets are freely available, and most training runs take only a few minutes. AI use statement We used AI tools to: • Conduct literature reviews and locate relevant papers. Every cited paper has been read by at least one co-author. • Write and maintain the code for this paper, including the ablation scripts. All code has been manually reviewed by two co-authors, and large parts were written by hand. • Draft and edit LaTeX, including the theory section and the proofs. We reviewed the text and rewrote most of it. We take full responsibility for the claims made in this paper. References Ahmed et al. (1974) N. Ahmed, T. Natarajan, and K. R. Rao Discrete cosine transform. IEEE transactions on Computers 100 (1), p. 90–93. Cited by: §B.3, §2.2. Alexandrov et al. (2020) A. Alexandrov, K. Benidis, M. Bohlke-Schneider, V. Flunkert, J. Gasthaus, T. Januschowski, D. C. Maddix, S. Rangapuram, D. Salinas, J. Schulz, et al. GluonTS: probabilistic and neural time series modeling in Python. Journal of Machine Learning Research 21 (116), p. 1–6. Cited by: §C.1, §C.2, §3. Baraniuk and Jones (1995) R. G. Baraniuk and D. L. Jones Unitary equivalence: a new twist on signal processing. IEEE Transactions on Signal Processing 43 (10), p. 2269–2282. Cited by: Appendix A, §2.2. Chan (1989) R. H. Chan Circulant preconditioners for Hermitian Toeplitz systems. SIAM Journal on Matrix Analysis and Applications 10 (4), p. 542–550. Cited by: Appendix A. Chen and Sun (2026) R. Chen and F. Sun DMamba: decomposition-enhanced mamba for time series forecasting. arXiv preprint arXiv:2602.09081. Cited by: Table 14, Table 15. Crabbé et al. (2024) J. Crabbé, N. Huynh, J. Stanczuk, and M. van der Schaar Time series diffusion in the frequency domain. arXiv preprint arXiv:2402.05933. Cited by: Appendix A, §1. Gabor (1946) D. Gabor Theory of communication. Journal of the Institution of Electrical Engineers – Part I 93 (26), p. 429–441. Cited by: Appendix A, §1, §2.2. Gray (2006) R. M. Gray Toeplitz and circulant matrices: a review. Foundations and Trends in Communications and Information Theory 2 (3), p. 155–239. Cited by: Appendix A, §B.3, §B.3, §2.2. Grenander and Szegő (1958) U. Grenander and G. Szegő Toeplitz forms and their applications. University of California Press. Cited by: Appendix A. Gu et al. (2022) A. Gu, K. Goel, and C. Ré Efficiently modeling long sequences with structured state spaces. In International Conference on Learning Representations, Cited by: Appendix A, §1. Gupta et al. (2018) V. Gupta, T. Koren, and Y. Singer Shampoo: preconditioned stochastic tensor optimization. In International Conference on Machine Learning, Cited by: Appendix A. Jolliffe (2005) I. Jolliffe Principal component analysis. Encyclopedia of statistics in behavioral science. Cited by: §2.2. Jordan et al. (2024) K. Jordan, Y. Jin, V. Boza, J. You, F. Cesista, L. Newhouse, and J. Bernstein Muon: an optimizer for hidden layers in neural networks. Note: https://kellerjordan.github.io/posts/muon/ Cited by: Appendix A. Kac et al. (1953) M. Kac, W. L. Murdock, and G. Szegő On the eigenvalues of certain Hermitian forms. Journal of Rational Mechanics and Analysis 2, p. 767–800. Cited by: Appendix A. Karhunen (1947) K. Karhunen Über lineare Methoden in der Wahrscheinlichkeitsrechnung. Ann. Acad. Sci. Fennicae. Ser. A. I. 37, p. 1–79. Cited by: Appendix A, §1, §1. Kingma and Ba (2015) D. P. Kingma and J. Ba Adam: a method for stochastic optimization. In International Conference on Learning Representations, Cited by: Appendix A, §2.4. Kollovieh et al. (2025) M. Kollovieh, M. Lienen, D. Lüdke, L. Schwinn, and S. Günnemann Flow matching with Gaussian process priors for probabilistic time series forecasting. In International Conference on Learning Representations, Cited by: Appendix A, §C.1, §C.2, Table 14, Table 15, §1, §2.3, §2, §3. Lipman et al. (2023) Y. Lipman, R. T. Q. Chen, H. Ben-Hamu, M. Nickel, and M. Le Flow matching for generative modeling. In International Conference on Learning Representations, Cited by: Appendix A, §1, §2.3. Loève (1978) M. Loève Probability theory I. 4 edition, Springer. Cited by: Appendix A. Martens and Grosse (2015) J. Martens and R. Grosse Optimizing neural networks with Kronecker-factored approximate curvature. In International Conference on Machine Learning, Cited by: §2.4. Mercer (1909) J. Mercer Functions of positive and negative type, and their connection with the theory of integral equations. Philosophical Transactions of the Royal Society of London. Series A 209, p. 415–446. Cited by: Appendix A, §1, §1. Nesterov (2004) Y. Nesterov Introductory lectures on convex optimization: a basic course. Springer. Cited by: §B.1, §1, §2.4, Proposition 1. Nie et al. (2022) Y. Nie, N. H. Nguyen, P. Sinthong, and J. Kalagnanam A time series is worth 64 words: long-term forecasting with transformers. arXiv preprint arXiv:2211.14730. Cited by: Table 14, Table 15. Ning et al. (2025) M. Ning, M. Li, J. Su, H. Jia, L. Liu, M. Beneš, W. Chen, A. A. Salah, and I. O. Ertugrul DCTdiff: intriguing properties of image generative modeling in the DCT space. In International Conference on Machine Learning, Cited by: Appendix A, §1. Phillips et al. (2022) A. Phillips, T. Seror, M. Hutchinson, V. De Bortoli, A. Doucet, and E. Mathieu Spectral diffusion processes. In NeurIPS 2022 Workshop on Score-Based Methods, Cited by: Appendix A, §1. Rasmussen and Williams (2006) C. E. Rasmussen and C. K. I. Williams Gaussian processes for machine learning. MIT Press. Cited by: Appendix A, §B.3, §2.2. Rasul et al. (2021) K. Rasul, C. Seward, I. Schuster, and R. Vollgraf Autoregressive denoising diffusion models for multivariate probabilistic time series forecasting. In International Conference on Machine Learning, Cited by: Appendix A. Tashiro et al. (2021) Y. Tashiro, J. Song, Y. Song, and S. Ermon CSDI: conditional score-based diffusion models for probabilistic time series imputation. In Advances in Neural Information Processing Systems, Cited by: Appendix A. Tong et al. (2024) A. Tong, K. Fatras, N. Malkin, G. Huguet, Y. Zhang, J. Rector-Brooks, G. Wolf, and Y. Bengio Improving and generalizing flow-based generative models with minibatch optimal transport. In International Conference on Artificial Intelligence and Statistics, Cited by: Appendix A, §1, §2.3. van der Sluis (1969) A. van der Sluis Condition numbers and equilibration of matrices. Numerische Mathematik 14, p. 14–23. Cited by: §2.4. Vyas et al. (2025) N. Vyas, D. Morwani, R. Zhao, M. Kwun, I. Shapira, D. Brandfonbrener, L. Janson, and S. Kakade SOAP: improving and stabilizing Shampoo using Adam. In International Conference on Learning Representations, Cited by: Appendix A. Whittle (1951) P. Whittle Hypothesis testing in time series analysis. Almqvist & Wiksell. Cited by: Appendix A. Xie et al. (2025) S. Xie, M. A. Mohamadi, and Z. Li Adam exploits ℓ∞ _∞-geometry of loss landscape via coordinate-wise adaptivity. In International Conference on Learning Representations, Cited by: Appendix A, §2.4. Zeng et al. (2023) A. Zeng, M. Chen, L. Zhang, and Q. Xu Are transformers effective for time series forecasting?. In Proceedings of the AAAI conference on artificial intelligence, Vol. 37, p. 11121–11128. Cited by: Table 14, Table 15. Zhang et al. (2025) T. Zhang, L. Maes, A. Milligan, A. Jolicoeur-Martineau, I. Mitliagkas, D. Scieur, S. Lacoste-Julien, and C. Guille-Escuret Understanding Adam requires better rotation dependent assumptions. In Advances in Neural Information Processing Systems, Cited by: Appendix A, §B.1, §1, §2.4, Proposition 1. Appendix A Related work Probabilistic generative forecasting. Denoising diffusion is a standard sampler for probabilistic time series: TimeGrad (Rasul et al., 2021) forecasts autoregressively with a diffusion decoder, and CSDI (Tashiro et al., 2021) performs conditional score matching for imputation and forecasting in observation coordinates. Conditional flow matching (Lipman et al., 2023; Tong et al., 2024) replaces the reverse SDE with regression onto a known interpolant velocity. TSFlow (Kollovieh et al., 2025) applies that recipe to forecasting with a Gaussian-process prior on the forecast window, Euclidean CFM in observation coordinates, and a bidirectional S4 backbone (Gu et al., 2022). We keep that prior and interpolant, and move the flow to a linear latent of the same window so that an MLP can replace the sequential backbone. Latent-space generative models. Several recent generative models run the denoising process in an orthonormal latent rather than in observation coordinates. DCTdiff (Ning et al., 2025) trains image diffusion in block DCT coordinates and reports gains in sample quality and training cost over pixel-space diffusion. Crabbé et al. (Crabbé et al., 2024) dualise time-series diffusion to the Fourier domain under mirrored Brownian motion, so that frequency-localised series are scored where they concentrate. The DCT is near Karhunen–Loève for smooth stationary Grams on a finite interval, and the DFT diagonalises every circulant covariance; Whittle and Szegő treat it as an approximate KL basis for stationary Toeplitz Grams (Gabor, 1946; Baraniuk and Jones, 1995; Whittle, 1951; Kac et al., 1953; Grenander and Szegő, 1958; Gray, 2006; Chan, 1989). Phillips et al. (Phillips et al., 2022) go further toward a kernel eigenbasis: they expand functional data in truncated Mercer coefficients and train a score model on the retained modes, with the kernel taken as the empirical covariance or a stationary alternative such as RBF or Matérn. Mercer’s theorem supplies that eigenbasis for a continuous positive-definite kernel (Mercer, 1909), and the Karhunen–Loève theorem yields uncorrelated coefficients for a second-order process, independent under a GP (Karhunen, 1947; Loève, 1978; Rasmussen and Williams, 2006). Although data-matched GP priors are natural for structured time series, existing approaches do not examine how such priors interact with latent projections under flow matching. Our contribution is to study this interaction: fixing a forecasting GP prior, we compare five full-order codecs on one MLP and one CFM loss. Rotation equivariance of optimisers. Euclidean SGD is equivariant to orthogonal reparameterisations of the weights: conjugating the network by a rotation conjugates the iterates and leaves the loss trajectory unchanged. Adam (Kingma and Ba, 2015) is not: its per-coordinate second-moment rescaling depends on the parameter basis, so rotating the coordinates changes the optimisation path (Zhang et al., 2025; Xie et al., 2025). That gap has motivated rotation-aware and approximately rotation-equivariant optimisers, including Kronecker-factored preconditioners such as Shampoo (Gupta et al., 2018), Shampoo–Adam hybrids such as SOAP (Vyas et al., 2025), and orthogonalised updates such as Muon (Jordan et al., 2024). Proposition 1 extends to any optimiser whose update commutes with W↦WQW WQ, which includes Muon and Shampoo without grafting; we test SGD only. Appendix B Theoretical Proofs B.1 Proof of Proposition 1 (SGD Invariance) Fix a square orthogonal codec matrix Φ∈ℝL×L ^L× L acting on the ambient observation vector ∈ℝLy ^L that enters the first linear layer. (For the concatenated joint window L=Lp+LfL=L_p+L_f, this is precisely Φ ; when context and horizon are encoded independently, one uses a block-diagonal isometry Q=blkdiag(Φp,Φf)Q=blkdiag( _p, _f), but the algebraic conjugation remains identical). Let 0=[p,0]z_0=[y_p,y_0] and 1=[p,f]z_1=[y_p,y_f] denote the ambient source and target joint states, respectively, yielding latent endpoints 0=Φ⊤0c_0= z_0 and 1=Φ⊤1c_1= z_1. The target displacement in coefficient space is =1−0=Φ⊤(1−0)u=c_1-c_0= (z_1-z_0). The conditional flow matching probability path is the image of the ambient observation-space line under Φ⊤ . Because Φ is an orthogonal isometry (ΦΦ⊤=IL =I_L), it preserves Euclidean norms: ‖2=‖Φ⊤(1−0)‖2=‖1−0‖2.\|u\|_2=\| (z_1-z_0)\|_2=\|z_1-z_0\|_2. (17) Consequently, the squared conditional flow matching regression objective equation 10 is isometric across all orthogonal coordinate systems. The hypothesis of Proposition 1 holds for any velocity field whose first and final linear projections bracket an entrywise hidden mapping whose internal parameters are independent of Φ . The implemented MLP architecture satisfies this condition (Section C.3). Denoting by ψ(⋅,t)ψ(·,t) the bracketed entrywise nonlinear operation (comprising hidden layers and AdaLN time modulation), the network computing the velocity field θv_θ has the form: =W1+1,~=ψ(,t,θψ),θ(t,)=W2~+2,h=W_1c+b_1, h=ψ(h,\,t;\, _ψ), _θ(t,c)=W_2 h+b_2, (18) where ψ(⋅,t)ψ(·,t) acts entrywise on its first argument, with parameters governed exclusively by the flow time t. Consider an Identity network parameterised by θI=W1,1,W2,2,θψ _I=\W_1,b_1,W_2,b_2, _ψ\ and define the conjugated parameter set θΦ=W1Φ,1Φ,W2Φ,2Φ,θψΦ _ =\W_1 ,b_1 ,W_2 ,b_2 , _ψ \ of the codec-Φ network by W1Φ=W1Φ,W2Φ=Φ⊤W2,2Φ=Φ⊤2,1Φ=1,θψΦ=θψ.W_1 =W_1 , W_2 = W_2, _2 = b_2, _1 =b_1, _ψ = _ψ. (19) Evaluating the conjugated velocity field on the rotated input =Φ⊤c= y yields: θΦ(t,Φ⊤) _ _ (t, y) =W2Φψ(W1ΦΦ⊤+1Φ,t,θψΦ)+2Φ =W_2 \,ψ (W_1 y+b_1 ,\,t;\, _ψ )+b_2 =Φ⊤W2ψ(W1ΦΦ⊤+1,t,θψ)+Φ⊤2 = W_2\,ψ (W_1 y+b_1,\,t;\, _ψ )+ b_2 =Φ⊤(W2ψ(W1+1,t,θψ)+2) = (W_2\,ψ(W_1y+b_1,\,t;\, _ψ)+b_2 ) =Φ⊤θI(t,). = v_ _I(t,y). (20) Subtracting the target displacement =Φ⊤(1−0)u= (z_1-z_0) confirms that the residual vectors are related by an exact orthogonal rotation: θΦ(t,Φ⊤)−=Φ⊤(θI(t,)−(1−0)).v_ _ (t, y)-u= (v_ _I(t,y)-(z_1-z_0) ). (21) Because the Euclidean norm is rotation-invariant, the instantaneous conditional flow matching loss is preserved identically at every point along the probability path: ‖θΦ−‖22=‖θI−(1−0)‖22.\|v_ _ -u\|_2^2=\|v_ _I-(z_1-z_0)\|_2^2. (22) Furthermore, the velocity field intertwines the decoded trajectories identically: if ˙=θΦ(t,) c=v_ _ (t,c) and =Φy= , then by the chain rule, ˙=Φ˙=ΦθΦ(t,Φ⊤)=θI(t,). y= c= _ _ (t, y)=v_ _I(t,y). (23) The equality of squared residuals establishes that the loss landscapes satisfy ℒΦ(θΦ)=ℒI(θI)L_ ( _ )=L_I( _I). Moreover, the parameter transformation equation 19 is an exact linear isometry under the standard Frobenius/Euclidean inner product on parameter space: ‖W1Φ‖F2+‖W2Φ‖F2+‖2Φ‖22+‖1Φ‖22=‖W1‖F2+‖W2‖F2+‖2‖22+‖1‖22.\|W_1 \|_F^2+\|W_2 \|_F^2+\|b_2 \|_2^2+\|b_1 \|_2^2=\|W_1\|_F^2+\|W_2\|_F^2+\|b_2\|_2^2+\|b_1\|_2^2. (24) Under rotation-invariant initialisation schemes, standard SGD and Euclidean momentum methods (such as heavy-ball momentum and Nesterov accelerated gradient) generate updates that are strictly equivariant to orthogonal conjugation (Nesterov, 2004; Zhang et al., 2025). That is, the gradient step on the conjugated parameters matches the conjugation of the Euclidean gradient step on the original parameters. Consequently, given conjugated initialisations, the training loss trajectories and decoded probability flows of Mercer, Fourier, DCT, PCA, and the unrotated Raw baseline coincide identically at every training iteration under Euclidean SGD. B.2 SGD Control on Real Benchmarks To probe whether the rotation-equivariance of SGD remains visible under our non-rotation-invariant initialisations and matched training budget, we re-train the five codecs with Nesterov SGD (momentum 0.90.9, learning rate 10−310^-3) on ETTh1 and ETTh2, leaving every other hyperparameter identical to the AdamW runs in the main text. The SGD learning rate is untuned; the comparison is the matched-budget control that Proposition 1 predicts. Table 3 reports test CRPS over seeds 4242–4444 at 2020 NFE under the shared joint window (L=384L=384) and OU length scales. All five codecs land within one standard deviation of each other on both datasets, in agreement with the equivariance claim. The shared SGD level is uniformly worse than the AdamW level on these two datasets (e.g. 0.22550.2255 vs. 0.1880.188 for Mercer on ETTh1, about 20%20\% above AdamW). On synthetic OU paths where the analytical prior matches the data-generating law exactly (Section F.1), SGD performs on par with AdamW (0.60330.6033 vs. 0.60630.6063), indicating that Adam’s coordinate-wise preconditioning primarily provides an optimisation advantage on complex, non-stationary data. Table 3: SGD ablation on ETTh1 and ETTh2 (Test CRPS, ↓ ) for Post-AdaLN Nesterov SGD (momentum 0.90.9, learning rate 10−310^-3). Evaluated at 2020 NFE under the identical joint window (L=384L=384) and OU length scales as the AdamW runs. Codec ETTh1 ETTh2 Mercer (Ours) 0.2255±0.00070.2255± 0.0007 0.145±0.0020.145± 0.002 Fourier 0.2240±0.00200.2240± 0.0020 0.145±0.0030.145± 0.003 DCT 0.2260±0.00100.2260± 0.0010 0.145±0.0020.145± 0.002 PCA 0.2248±0.00040.2248± 0.0004 0.145±0.0020.145± 0.002 Raw 0.2240±0.00100.2240± 0.0010 0.145±0.0030.145± 0.003 B.3 Proof of Proposition 2 (Jacobi Preconditioning) Let A=[⊤]≻0A=E[xx ] 0 under the centred prior on the ambient grid, and set D=diag(A)−1/2D=diag(A)^-1/2. Under the standing model that Adam’s input-side rescaling acts as this Jacobi preconditioner, the change of coordinates ↦Dx Dx leaves the residual Gram G=DAD,Gij=AijAiiAjj.G=DAD, G_ij= A_ij A_iA_j. Its condition number is κ(G)=λmax(G)/λmin(G)κ(G)= _ (G)/ _ (G), as in equation 14. The spectrum of A is invariant under any orthogonal codec ↦Q⊤x Q x, since Q⊤AQ AQ is similar to A. The diagonal of A is not. The uncentred conditional interpolant used in training is outside the claim (cf. the measurements in §3). The proposition’s summary is given in Table 4, built by evaluating the closed form of A for each codec, then applying the Jacobi formula D=diag(A)−1/2D=diag(A)^-1/2. The calculations below verify each row. Table 4: Theoretical first-layer Gram matrix A=[⊤]A=E[xx ] and corresponding Jacobi condition number κ(DAD)κ(DAD) under the centred Ornstein–Uhlenbeck prior ∼(,K)y (0,K) on a window of length n=Ln=L, with regularizing jitter ε>0 >0. Codec First-layer Gram Matrix A Jacobi Condition Number Bound κ(DAD)κ(DAD) Identity (Raw) Tn(f)+εInT_n(f)+ I_n κ(A)κ(A), increases with length scale ℓ Mercer (Ours) Λ=diag(λ1,…,λn) =diag( _1,…, _n) 11 (exact decorrelation) Fourier ΦF⊤Tn(f)ΦF+εIn _F T_n(f)\, _F+ I_n 1+(1/minωf(ω))1+O (1/ _ωf(ω) ) DCT ΦD⊤Tn(f)ΦD+εIn _D T_n(f)\, _D+ I_n 1+(1/minωf(ω))1+O (1/ _ωf(ω) ) PCA Empirical data covariance Data-dependent (not aligned with prior kernel) Identity. A stationary kernel on a regular grid makes Cov()Cov(x) Toeplitz, hence constant on the diagonal: Ajj=a:=k(0)+εA_j=a:=k(0)+ for every j. Then D=a−1/2ID=a^-1/2I and G=A/aG=A/a, so κ(G)=κ(A)κ(G)=κ(A). Mercer. With ∼(,K)y (0,K), K=ΦΛΦ⊤K= , and Φ square orthogonal, =Φ⊤⇒A=[⊤]=Φ⊤KΦ=Λc= y A=E[cc ]= K = by Mercer’s theorem and the Karhunen–Loève expansion (Rasmussen and Williams, 2006). Jitter gives λi>0 _i>0, so D=Λ−1/2D= ^-1/2 and G=IG=I, hence κ(G)=1κ(G)=1. Fourier. On a window of length n, K=Tn(f)K=T_n(f) is the symmetric Toeplitz covariance generated by the spectral density f. Let CnC_n be the optimal Strang circulant approximation to Tn(f)T_n(f), which is exactly diagonalised by the discrete Fourier basis: ΦF⊤CnΦF=Δ:=diag(f(2πk/n)) _F C_n _F= :=diag(f(2π k/n)) with Δ≻0 0 (Gray, 2006). The latent second-moment matrix is A=ΦF⊤Tn(f)ΦFA= _F T_n(f) _F, whose diagonal asymptotically satisfies Akk=f(2πk/n)+O(n−1)A_k=f(2π k/n)+O(n^-1). Under the Jacobi preconditioner D=diag(A)−1/2D=diag(A)^-1/2, the preconditioned Gram evaluates to: DAD−I=DΦF⊤(Tn(f)−Cn)ΦFD+O(n−1).DAD-I=D _F (T_n(f)-C_n) _FD+O(n^-1). Because Toeplitz covariance does not wrap periodically, the boundary discrepancy Tn−CnT_n-C_n is a dense triangular perturbation driven by the edge jump k(0)−k(n−1)=1−e−(n−1)/ℓk(0)-k(n-1)=1-e^-(n-1)/ . In the Fourier domain, this step discontinuity generates harmonic spectral leakage (the Gibbs phenomenon). Because the Jacobi weights scale as Dkk≈f(2πk/n)−1/2D_k≈ f(2π k/n)^-1/2, the off-diagonal entries scale inversely with high-frequency spectral density: (DAD)jk∼⟨ϕj,(Tn−Cn)ϕk⟩f(ωj)f(ωk).(DAD)_jk _j,(T_n-C_n) _k f( _j)f( _k). For an OU process, f(ω)∝ℓ−1ℓ−2+ω2f(ω) ^-1 ^-2+ω^2. When ℓ is large, the spectral mass concentrates at low frequencies, so minkf _kf shrinks rapidly. Dividing by the vanishing high-frequency spectral density severely amplifies the off-diagonal elements of DADDAD, inflating κ(DAD)κ(DAD). The empirical numbers in Table 5 precisely follow this dependency: κ(DAD)κ(DAD) scales monotonically with ℓ , growing from 13.413.4 on Traffic (ℓ=7 =7) to 268.5268.5 on ETTh1 (ℓ=336 =336). DCT. The DCT-I basis ΦD _D diagonalises the symmetric circulant CnsymC_n^sym that approximates Tn(f)T_n(f) under mirror-reflective boundary conditions (Ahmed et al., 1974; Gray, 2006). Crucially, even reflection across the grid boundaries eliminates the periodic wrap-around step discontinuity k(0)−k(n−1)k(0)-k(n-1) that plagues the Fourier basis. Consequently, the boundary discrepancy Tn(f)−CnsymT_n(f)-C_n^sym produces negligible spectral leakage, leaving the off-diagonal entries of ΦD⊤TnΦD _D T_n _D tightly suppressed across all frequencies. Under Jacobi scaling, DAD≈I+EsymDAD≈ I+E_sym, where the perturbation ‖Esym‖2\|E_sym\|_2 remains uniformly small and stable with respect to the length scale ℓ . As reported in Table 5, DCT maintains a near-optimal condition number κ(DAD)∈[1.90,2.12]κ(DAD)∈[1.90,2.12] across all benchmarks. PCA. PCA’s basis is the eigenbasis of the empirical training covariance, which is data-dependent. On the centred prior population of the proposition, PCA’s basis is generally not aligned with K, and κ(DAD)κ(DAD) is whatever that misalignment produces; Table 5 reports 8.158.15–46.9846.98. On its own empirical-covariance population A=diag(λiPCA)A=diag(λ^PCA_i) the codec gives κ(DAD)=1κ(DAD)=1 trivially; this is a separate statement and is not the population the proposition addresses. Empirical check. The population the identity is proved against is the centred prior y∼(0,K)y (0,K) on the joint grid; the training input is not that population. Table 5 reports κ(DAD)κ(DAD) on it: Mercer hits 11 on every dataset, DCT sits at 1.91.9–2.12.1 (reflecting boundary reflection stability), Fourier in the tens to hundreds (growing with ℓ ), PCA at 88–4747, and Raw from 1.96×1021.96× 10^2 on Traffic to 1.72×1051.72× 10^5 on ETTh1. Table 5: Empirical Jacobi condition number κ(DAD)κ(DAD) under the centred prior y∼(0,K)y (0,K) on the joint grid (joint window 192192/192192; last channel; scaler fit on the first 70%70\%). This is the population the proposition’s identity is proved against: Mercer hits κ(DAD)=1κ(DAD)=1 on every dataset, with jitter only, by construction. The training input is not this population; Table 12 reports the same quantity under the conditional draw 0=Φ⊤[p;0]c_0= [y_p;y_0]. Dataset Raw Fourier Mercer DCT PCA ETTh1 1.72×1051.72× 10^5 268.50268.50 1.001.00 1.921.92 8.158.15 ETTh2 8.16×1048.16× 10^4 189.71189.71 1.001.00 2.072.07 46.9846.98 Weather 2.80×1042.80× 10^4 119.74119.74 1.001.00 2.122.12 13.4613.46 Electricity 2.23×1032.23× 10^3 41.6341.63 1.001.00 2.082.08 43.9343.93 Traffic 195.79195.79 13.4413.44 1.001.00 1.901.90 44.6844.68 Appendix C Experimental Setup and Implementation Details C.1 Evaluation Metrics To evaluate forecast distribution accuracy and point calibration, we draw S=100S=100 sample trajectories ^f(s)s=1S\ y_f^(s)\_s=1^S per test window in the observation space (after inverse scaling). Continuous Ranked Probability Score (CRPS). Following standard probabilistic forecasting benchmarks (Alexandrov et al., 2020; Kollovieh et al., 2025), we report the normalised CRPS, defined as the mean weighted quantile loss across nine standard evaluation quantiles =0.1,0.2,…,0.9Q=\0.1,0.2,…,0.9\: CRPS=∑q∈∑j=0Lf−1ρq(yf,j−y^f,j(q))||∑j=0Lf−1|yf,j|,CRPS= _q _j=0^L_f-1 _q (y_f,j- y_f,j^(q) )|Q| _j=0^L_f-1|y_f,j|, (25) where y^f,j(q) y_f,j^(q) denotes the empirical q-quantile of the S sampled forecasts at horizon step j, and ρq(e)=2max(qe,(q−1)e) _q(e)=2 (qe,(q-1)e) is the asymmetric pinball loss. Mean Absolute Error (MAE). For deterministic point forecast evaluation, we take the empirical median across generated trajectory samples, ^fmed=median(^f(s)s=1S) y_f^med=median(\ y_f^(s)\_s=1^S), and evaluate the unnormalized absolute deviation: MAE=1Lf∑j=0Lf−1|yf,j−y^f,jmed|.MAE= 1L_f _j=0^L_f-1 |y_f,j- y_f,j^med |. (26) Both metrics are averaged across various rolling test horizons. C.2 Experimental Protocol All time-series channels are treated as independent univariate sequences. The observation-space prior is the Gaussian Process defined in Section 2.1, governed by the stationary Ornstein–Uhlenbeck (OU) covariance kernel with length scale ℓ selected via validation-likelihood maximisation (Table 8) and shared across all five MLP codecs. The context past and forecast horizon lengths are fixed at Lp=Lf=192L_p=L_f=192, yielding a total sequence window L=Lp+Lf=384L=L_p+L_f=384. Source prior samples and target data trajectories are paired independently during training; minibatch optimal transport is not required due to the structured probability paths induced by the GP prior. Joint Window Conditioning. The complete joint sequence of length L=384L=384 is mapped into the latent coefficient space at once via the square orthonormal matrix Φ∈ℝL×L ^L× L. The lightweight MLP described in Section C.3 maps the joint coefficient vector c∈ℝLc ^L to a hidden dimension of width 6464 (a 22-layer MLP with 11 hidden layer, ≈55k≈ 55k parameters), applies adaptive flow-time conditioning, and projects back to predict the velocity field. Parameter counts across Mercer, Fourier, Raw, DCT, and PCA are strictly matched. TSFlow Baseline. The reference flow matching baseline is TSFlow (Kollovieh et al., 2025), evaluated on the identical context-and-horizon window (Lp=Lf=192L_p=L_f=192) and shared length scale ℓ . TSFlow operates in ambient observation space and utilises an expressive bidirectional S4 backbone comprising 186,833186,833 parameters. Evaluation Protocol. Predicted distributions are evaluated using the normalised Continuous Ranked Probability Score (CRPS) across nine evaluation quantiles Q=0.1,…,0.9Q=\0.1,…,0.9\ (Alexandrov et al., 2020), normalised by the sum of absolute ground-truth horizon values scaled by |Q||Q|, i.e., |Q|∑j=0Lf−1|yf,j||Q| _j=0^L_f-1|y_f,j|. Deterministic point forecasts are evaluated via the Mean Absolute Error (MAE) of the empirical median trajectory ^fmed y_f^med. All metrics are computed via 2020 Euler integration steps (NFE) with S=100S=100 sample trajectories per test window, averaged over random seeds 4242–4848. C.3 AdaLN Architectural Placement To ensure our network satisfies the conditions of Proposition 1, the implemented MLP applies adaptive time modulation (in the style of AdaLN) after the first linear layer. Given the concatenated input vector ∈ℝLc ^L (joint sequence coefficients) and the sinusoidal flow-time embedding t∈ℝdee_t ^d_e, the forward pass is defined as: =W1+1,~=GELU((1+γ(t))⊙+β(t)),θ=W2~+2,h=W_1c+b_1, h=GELU ((1+γ(e_t)) +β(e_t) ), _θ=W_2 h+b_2, (27) where the scale and shift vectors γ(t),β(t)∈ℝdhγ(e_t),β(e_t) ^d_h are generated via a zero-initialised linear projection of te_t. The bracketed operation ψ(,t)=GELU((1+γ(t))⊙+β(t))ψ(h,t)=GELU ((1+γ(e_t)) +β(e_t) ) (28) acts purely entrywise on the hidden activation vector h, with modulation parameters depending exclusively on flow time t. This strictly satisfies the structural condition required in the proof of Proposition 1 (Appendix B.1). C.4 Dataset Properties and Validation Likelihood Each benchmark time series is assigned the covariance kernel achieving the lowest Gaussian negative log-likelihood on its validation slice. A standard scaler is fit on the first 70%70\% of the series. The validation score is computed as the mean negative log density across the subsequent 15%15\% under a zero-mean Gaussian process. For the stationary Ornstein–Uhlenbeck (OU) and squared exponential (RBF) kernels, the continuous length scale is searched over ℓ∈[0.3,2000] ∈[0.3,2000]. For the periodic kernel, the length scale is searched over ℓ∈[0.02,8] ∈[0.02,8] alongside an integer period p∈[4,min(400,Nval/2)]p∈[4, (400,N_val/2)]. Every candidate kernel includes a regularizing diagonal jitter ε∈[10−5,1.0] ∈[10^-5,1.0]. Hyperparameters are optimized using 400400 Tree-structured Parzen Estimator (TPE) trials. As an intuitive baseline, isotropic white noise on a series scaled to unit variance scores approximately 1.421.42. Table 6 compares the validation negative log-likelihood across candidate kernels. Traffic exhibits dominant weekly periodicity (p=168p=168). ETTh2 and Electricity favor the smooth squared exponential (RBF) kernel, whereas ETTh1 and Weather are optimal under the rougher Ornstein–Uhlenbeck (OU) kernel. Table 7 details the optimal continuous hyperparameters found via TPE for each candidate kernel. In our primary benchmark experiments, we adopt standard domain-aligned integer length scales ℓ for the baseline OU prior (matching characteristic temporal cycles close to the continuous validation optima), along with a runtime diagonal jitter of ε=10−4 =10^-4 shared across all models; these parameters are summarized in Table 8. Table 6: Mean negative log density of candidate Gaussian process priors on the validation slice. Lower is better (↓ ). The periodic entry for Traffic uses period p=168p=168, length scale ℓ=0.1 =0.1, and jitter ε=0.05 =0.05; all other entries represent the optimal fit across 400400 TPE trials. Bold indicates the best score in each row. Dataset OU RBF Periodic ETTh1 −1.231-1.231 −0.970-0.970 0.0620.062 ETTh2 −0.836-0.836 −1.184-1.184 0.6080.608 Traffic 0.4650.465 0.2260.226 0.0350.035 Electricity 0.2280.228 0.2090.209 0.9690.969 Weather −0.587-0.587 −0.372-0.372 1.2031.203 Table 7: Continuous optimal parameters of the validation likelihood fits from Table 6. Length scale ℓ , period p (periodic kernel only), and regularizing diagonal jitter ε . The jitter search interval is ε∈[10−5,1.0] ∈[10^-5,1.0]; the OU fits for ETTh2, Traffic, Electricity, and Weather saturate at the lower bound 1.0×10−51.0× 10^-5. Dataset Kernel ℓ Period p Jitter ε ETTh1 OU 479479 — 4.4×10−44.4× 10^-4 ETTh1 RBF 7.417.41 — 2.8×10−32.8× 10^-3 ETTh1 Periodic 1.981.98 276276 0.0620.062 ETTh2 OU 183183 — 1.0×10−51.0× 10^-5 ETTh2 RBF 4.584.58 — 7.2×10−47.2× 10^-4 ETTh2 Periodic 2.772.77 2424 0.230.23 Traffic OU 13.113.1 — 1.1×10−51.1× 10^-5 Traffic RBF 3.703.70 — 0.0250.025 Traffic Periodic 0.100.10 168168 0.0500.050 Electricity OU 20.720.7 — 1.1×10−51.1× 10^-5 Electricity RBF 3.523.52 — 0.0280.028 Electricity Periodic 0.460.46 168168 0.530.53 Weather OU 111111 — 1.1×10−51.1× 10^-5 Weather RBF 5.575.57 — 0.0100.010 Weather Periodic 0.670.67 288288 0.640.64 Table 8: Baseline OU length scale ℓ and runtime diagonal jitter ε used in primary benchmark runs, shared across Mercer, Fourier, DCT, PCA, and Raw. Dataset Length Scale ℓ Runtime Jitter ε ETTh1 336336 10−410^-4 ETTh2 192192 10−410^-4 Weather 9696 10−410^-4 Electricity 2424 10−410^-4 Traffic 77 10−410^-4 First-Layer Conditioning under Traffic Covariance Kernels. When evaluating linear codecs on Traffic, departing from the stationary Ornstein–Uhlenbeck (OU) assumption breaks the near-Toeplitz covariance structure that DCT diagonalizes. Table 9 reports the conditional Jacobi condition number κ(DAD)κ(DAD) evaluated during training under both the baseline OU prior (ℓ=7 =7, ε=10−4 =10^-4) and the optimal weekly periodic prior (p=168p=168, ℓ=0.1 =0.1, ε=0.05 =0.05). While DCT provides a competitive condition number under the OU prior (κ(DAD)=712κ(DAD)=712), its condition number surges to 2.45×1042.45× 10^4 under the periodic kernel, confirming that fixed harmonic bases cannot accommodate non-Toeplitz periodic correlations. In contrast, MercerFlow analytically derives its orthonormal basis directly from the kernel Gram matrix K, dynamically adapting to periodic covariance and maintaining a tightly bounded condition number (κ(DAD)=1.18×103κ(DAD)=1.18× 10^3). Table 9: Traffic first-layer Jacobi condition number κ(DAD)κ(DAD) under conditional prior samples c0=Φ⊤[yp,y0]c_0= [y_p,y_0]: stationary OU prior (ℓ=7 =7, ε=10−4 =10^-4) versus validation-optimal periodic prior (period p=168p=168, ℓ=0.1 =0.1, ε=0.05 =0.05). Evaluated on the joint window (Lp=192,Lf=192L_p=192,L_f=192, L=384L=384), last channel, random seed 4242, standard scaler fit on the first 70%70\%. Bold indicates the best (lowest) condition number in each row. Prior Mercer Fourier Raw DCT PCA OU 687687 891891 8.54×1038.54× 10^3 712712 477 Periodic 1.18×1031.18× 10^3 188188 1.62×1041.62× 10^4 2.45×1042.45× 10^4 23.023.0 Forecasting Performance Beyond Stationary Priors. To examine how prior–codec alignment translates to empirical forecasting quality, Table 10 evaluates all five linear codecs on Traffic under both the baseline OU prior and the validation-optimal weekly periodic prior at the default learning rate (η=10−3η=10^-3) and under learning rate stress (η=10−1η=10^-1). At the default learning rate (η=10−3η=10^-3), adopting the kernel-matched periodic prior improves CRPS across all five codecs, confirming that a data-aligned GP prior benefits flow matching globally. MercerFlow achieves the best overall forecast (0.08210.0821 CRPS). Crucially, when pushed to an aggressive learning rate (η=10−1η=10^-1), fixed harmonic codecs that assume stationary Toeplitz Gram matrices degrade severely: DCT CRPS deteriorates to 0.16410.1641, while Fourier becomes numerically unstable (CRPS explodes to 0.77730.7773). In contrast, MercerFlow preserves exact Karhunen–Loève decorrelation for the periodic kernel and maintains stable convergence (0.09510.0951 CRPS), demonstrating robust universality beyond stationary priors. Fourier Degradation Across Stationary and Periodic Regimes. Although Fourier captures periodic harmonic frequencies better than DCT on Traffic, it remains inferior to MercerFlow. On stationary processes with long temporal dependencies (such as ETTh1 and ETTh2), the finite-window boundary jump k(0)−k(L−1)k(0)-k(L-1) causes harmonic spectral leakage under Fourier coordinates. Because Adam’s coordinate preconditioning scales inversely with spectral density (Dkk∼f(ωk)−1/2D_k f( _k)^-1/2), dividing by the vanishing high-frequency spectral density severely inflates off-diagonal entries, driving Fourier’s condition number to 268.5268.5 on ETTh1 (Table 5) and 471471 under conditional training draws (Table 12), compared to 1.001.00 and 14.814.8 for MercerFlow. Under learning rate stress (η=10−1η=10^-1), these amplified high-frequency coordinates receive disproportionate gradient updates along leakage directions, making Fourier numerically unstable and causing seed divergence (0.77730.7773 CRPS in Table 10). By deriving its basis analytically from the exact finite-domain covariance, MercerFlow maintains exact Karhunen–Loève decorrelation without boundary leakage, ensuring stable convergence across learning rate regimes. Table 10: Traffic probabilistic forecasting performance (Test CRPS, 2020 NFE, ↓ ) under the stationary OU prior vs. the validation-optimal periodic prior (p=168p=168, ℓ=0.1 =0.1, ε=0.05 =0.05). Under the periodic prior, MercerFlow dynamically adapts and achieves top performance, whereas fixed harmonic codecs (DCT, Fourier) suffer from severe condition degradation and numerical instability under learning rate stress (η=10−1η=10^-1). Bold indicates the best score in each row. Prior Learning Rate η Mercer Fourier Raw DCT PCA OU 10−310^-3 0.0875±0.00090.0875± 0.0009 0.0871±0.00150.0871± 0.0015 0.0884±0.00030.0884± 0.0003 0.0878±0.00110.0878± 0.0011 0.0883±0.00180.0883± 0.0018 OU 10−110^-1 0.1065±0.00360.1065± 0.0036 0.1080±0.00860.1080± 0.0086 0.2462±0.07170.2462± 0.0717 0.1116±0.00260.1116± 0.0026 0.1065±0.00410.1065± 0.0041 Periodic 10−310^-3 0.0821±0.00060.0821± 0.0006 0.0827±0.00110.0827± 0.0011 0.0839±0.00160.0839± 0.0016 0.0840±0.00070.0840± 0.0007 0.0828±0.00070.0828± 0.0007 Periodic 10−110^-1 0.0951±0.01300.0951± 0.0130 0.7773±1.03500.7773± 1.0350 1.0198±0.10781.0198± 0.1078 0.1641±0.06880.1641± 0.0688 0.0855±0.00320.0855± 0.0032 C.5 Dataset Hardness and First-Layer Conditioning To quantify dataset difficulty independently of model architecture or codec choice, we define dataset “hardness” as the mean squared error (MSE) of the frozen analytical OU conditional mean evaluated on the held-out final 15%15\% slice (after standard scaling, averaged across all channels), using the optimal validation length scale ℓ from Table 8. Table 11 reports this residual alongside the validation negative log-likelihood. ETTh2 and Weather yield low residual errors of 0.330.33 and 0.440.44, respectively; we classify these as easy benchmarks where the analytical prior mean already captures the dominant temporal structure. In contrast, ETTh1, Electricity, and Traffic exhibit substantial residuals ranging from 1.001.00 to 1.461.46; we classify these as hard benchmarks where the velocity field must actively learn complex residual dynamics. Table 11: Dataset hardness from the frozen OU prior. Residual is the mean squared error of the conditional OU mean on the held-out last 15%15\%, after the training standard scaler, averaged over channels; length scale from Table 8. OU NLL is the validation score from Table 6; white noise on unit-scaled series scores about 1.421.42. Easy: residual <1/2<1/2; hard: residual ≥1≥ 1. ETTh2 Weather Electricity ETTh1 Traffic OU NLL −0.836-0.836 −0.587-0.587 0.2280.228 −1.231-1.231 0.4650.465 Residual 0.330.33 0.440.44 1.001.00 1.441.44 1.461.46 Empirical First-Layer Conditioning κ(DAD)κ(DAD). Table 12 reports the empirical Jacobi condition number κ(DAD)κ(DAD) evaluated on conditional prior samples 0=Φ⊤[p,0]c_0= [y_p,y_0] across both the training split and the held-out test split. On the training split, PCA achieves low condition numbers by construction, as its basis directly diagonalises the empirical training covariance. However, under test-set non-stationarity and distribution shift (notably on Electricity and Traffic), PCA’s condition number degrades sharply (surging from 115115 to 285285 on Electricity, and from 481481 to 12901290 on Traffic). Because MercerFlow derives its orthonormal basis analytically from the true prior kernel K, it exhibits consistent preconditioning across splits without empirical overfitting. Table 12: Empirical Jacobi condition number κ(DAD)κ(DAD) under the conditional prior draws 0=Φ⊤[p;0]c_0= [y_p;y_0] on the training split and on the held-out last 15%15\%, standard-scaled; D=diag(A)−1/2D=diag(A)^-1/2 with A the empirical second moment of c0c_0. Each cell gives train/test. Scaler fit on the first 70%70\%; the PCA basis is fit on data windows, not on prior draws. Raw values share a single exponent per row. Dataset Raw Fourier Mercer DCT PCA ETTh1 (2.08/3.25)(2.08/3.25)\,1e5 539539 / 471471 34.734.7 / 14.814.8 34.734.7 / 14.614.6 26.826.8 / 15.115.1 ETTh2 (3.45/6.21)(3.45/6.21)\,1e5 17701770 / 29702970 155155 / 248248 154154 / 260260 121121 / 198198 Weather (4.20/1.11)(4.20/1.11)\,1e4 295295 / 272272 13.513.5 / 4.54.5 14.314.3 / 4.64.6 17.617.6 / 15.915.9 Electricity (5.85/5.22)(5.85/5.22)\,1e3 234234 / 394394 60.260.2 / 171171 59.659.6 / 176176 115115 / 285285 Traffic (8.60/1.10)(8.60/1.10)\,1e3 922922 / 11401140 683683 / 925925 709709 / 936936 481481 / 12901290 Correlation Between Conditioning and Predictive Quality. To verify the direct connection between first-layer AdamW preconditioning and final forecasting quality, Table 13 reports the mean within-dataset Spearman rank correlation ρ between κ(DAD)κ(DAD) and test CRPS across the five MLP codecs. On the hard datasets, where the flow model must resolve significant residual dynamics beyond the prior mean, the rank correlation reaches ρ=0.79ρ=0.79 on held-out splits. On the easy datasets, the frozen prior already accounts for almost all predictive capacity, causing performance across codecs to saturate and making CRPS independent of first-layer conditioning (ρ=−0.08ρ=-0.08). Table 13: Mean across datasets of within-dataset Spearman ρ between κ(DAD)κ(DAD) under the conditional prior (Table 12) and CRPS, over the five MLP codecs; TSFlow excluded. Each per-dataset ρ ranks the five codecs on the two axes within that dataset; the column is the mean of those ρ’s across the listed datasets. n is the number of datasets averaged. κ(DAD)κ(DAD) Overall (n=5n=5) Hard (n=3n=3) Easy (n=2n=2) Test 0.440.44 0.790.79 −0.08-0.08 Train 0.380.38 0.680.68 −0.08-0.08 Appendix D Full Predictive Performance While the main text condenses benchmark evaluations into summary ranks and computational footprint (Section 3.1), this section provides comprehensive, uncompressed per-dataset results for both probabilistic and point forecasting across all five GluonTS benchmarks. D.1 Probabilistic Forecasting Performance (CRPS) Table 14 details the normalised Continuous Ranked Probability Score (CRPS, ×10−2× 10^-2) across all models, evaluated at 2020 Euler integration steps with 100100 sampled paths per forecast window. Parenthetical values indicate the model’s rank determined by pairwise two-sided t-tests on shared random seeds (p<0.2p<0.2) within the flow matching family. As shown, deterministic forecasters (DLinear, PatchTST, DMamba) exhibit severe CRPS degradation across every dataset, as they fail to capture multi-modal predictive distributions and uncertainty bounds. Among conditional flow matching models, MercerFlow consistently matches or outperforms all alternatives, securing rank 11 across all five benchmarks. Table 14: Comprehensive long-horizon probabilistic forecasting performance (Test CRPS, 10−210^-2, ↓ ). Top: deterministic baselines; bottom: conditional flow matching models. Parenthetical values denote paired two-sided t-test rank (p<0.2p<0.2) among flow matching models on shared random seeds (4242–4848). Bold: best; second best. Method ETTh1 ETTh2 Weather Electricity Traffic DLinear (Zeng et al., 2023) 23.06±0.0123.06± 0.01 16.02±0.0116.02± 0.01 6.26±0.026.26± 0.02 5.57±0.005.57± 0.00 10.75±0.0110.75± 0.01 PatchTST (Nie et al., 2022) 23.02±0.6123.02± 0.61 16.50±0.4716.50± 0.47 5.99±0.585.99± 0.58 6.70±0.206.70± 0.20 12.45±0.3012.45± 0.30 DMamba (Chen and Sun, 2026) 25.38±0.8925.38± 0.89 16.00±0.4116.00± 0.41 2.56±0.022.56± 0.02 5.48±0.105.48± 0.10 10.42±0.2410.42± 0.24 TSFlow (S4) (Kollovieh et al., 2025) 19.1±0.219.1± 0.2 (5) 15.4±0.315.4± 0.3 (6) 2.18±0.082.18± 0.08 (6) 7.4±0.37.4± 0.3 (6) 23.0±9.023.0± 9.0 (6) Raw 19.2±0.219.2± 0.2 (5) 13.1 ± 0.2 (1) 1.82±0.031.82± 0.03 (1) 4.23±0.064.23± 0.06 (5) 8.9±0.18.9± 0.1 (5) Fourier 19.0±0.319.0± 0.3 (1) 13.2±0.313.2± 0.3 (1) 1.82±0.031.82± 0.03 (1) 4.11±0.044.11± 0.04 (2) 8.70 ± 0.09 (1) DCT 18.9±0.318.9± 0.3 (1) 13.3±0.413.3± 0.4 (1) 1.82±0.021.82± 0.02 (1) 4.06 ± 0.08 (1) 8.8±0.18.8± 0.1 (3) PCA 18.9±0.118.9± 0.1 (1) 13.1 ± 0.4 (1) 1.81 ± 0.03 (1) 4.16±0.064.16± 0.06 (4) 8.8±0.18.8± 0.1 (2) MercerFlow (Ours) 18.8 ± 0.2 (1) 13.3±0.413.3± 0.4 (1) 1.82±0.031.82± 0.03 (1) 4.07±0.084.07± 0.08 (1) 8.7 ± 0.1 (1) D.2 Point Forecasting Performance (MAE) To assess point-forecasting capability, Table 15 reports the unnormalised Mean Absolute Error (MAE) evaluated on the empirical median of generated sample paths ^fmed=median(^f(s)s=1S) y_f^med=median(\ y_f^(s)\_s=1^S). Deterministic baselines often achieve strong MAE on smooth benchmarks (e.g., ETTh1 and ETTh2), but this comes at the expense of severe distributional degradation (Table 14). Conversely, MercerFlow delivers competitive point accuracy across all datasets while simultaneously preserving well-calibrated prediction intervals, maintaining top rank (1.001.00) among generative models. Table 15: Point forecasting performance (Test MAE, ↓ , Traffic in 10−310^-3). Parenthetical values denote paired two-sided t-test rank (p<0.2p<0.2) among flow matching models on shared random seeds (4242–4848). Bold: best; second best. Method ETTh1 ETTh2 Weather Electricity Traffic DLinear (Zeng et al., 2023) 2.12±0.012.12± 0.01 4.91±0.014.91± 0.01 27.56±0.0827.56± 0.08 187.2±0.2187.2± 0.2 4.00±0.014.00± 0.01 PatchTST (Nie et al., 2022) 2.11 ± 0.06 5.06±0.155.06± 0.15 26.38±2.5426.38± 2.54 225.4±6.8225.4± 6.8 5.00±0.015.00± 0.01 DMamba (Chen and Sun, 2026) 2.33±0.082.33± 0.08 4.91 ± 0.12 11.26±0.1011.26± 0.10 184.2±3.5184.2± 3.5 4.00±0.014.00± 0.01 TSFlow (S4) (Kollovieh et al., 2025) 2.25±0.032.25± 0.03 (1) 6.0±0.16.0± 0.1 (6) 11.6±0.311.6± 0.3 (6) 320±10320± 10 (6) 11.0±5.011.0± 5.0 (6) Raw 2.31±0.022.31± 0.02 (6) 5.13±0.095.13± 0.09 (1) 10.0 ± 0.2 (1) 181±3181± 3 (5) 4.04±0.084.04± 0.08 (3) Fourier 2.27±0.032.27± 0.03 (3) 5.1±0.15.1± 0.1 (1) 10.1±0.210.1± 0.2 (1) 176±2176± 2 (2) 3.95 ± 0.04 (1) DCT 2.25±0.042.25± 0.04 (1) 5.1±0.15.1± 0.1 (1) 10.0 ± 0.1 (1) 174 ± 3 (1) 4.02±0.054.02± 0.05 (3) PCA 2.26±0.022.26± 0.02 (1) 5.0 ± 0.1 (1) 10.0 ± 0.1 (1) 178±2178± 2 (4) 4.00±0.064.00± 0.06 (3) MercerFlow (Ours) 2.24±0.032.24± 0.03 (1) 5.1±0.15.1± 0.1 (1) 10.1±0.210.1± 0.2 (1) 174 ± 3 (1) 3.97±0.053.97± 0.05 (1) D.3 Sensitivity to Extended Forecasting Horizons To evaluate model resilience across extended horizons, we expand the forecast window from the default Lf=192L_f=192 to Lf∈336,720L_f∈\336,720\ while keeping the context length fixed at Lp=192L_p=192 (resulting in total sequence lengths of L=528L=528 and L=912L=912). Table 16 provides the full test CRPS and MAE comparisons across all five benchmarks. The results reveal distinct structural behaviors across datasets. On smooth series with long temporal correlation scales, such as ETTh1 (ℓ=336 =336), ambient observation-space modeling (Raw) remains effective. However, on complex, multi-frequency series such as Weather and Electricity at the longest horizon (Lf=720L_f=720), predicting directly in observation space incurs severe error accumulation and variance inflation (e.g., Raw MAE on Electricity degrades to 206.9±15.6206.9± 15.6, compared to 175.7±5.6175.7± 5.6 for DCT). In contrast, orthogonal latent codecs constrain error compounding and maintain calibration even as the window length quadruples, demonstrating that compact MLP architectures operating in decorrelated coordinates scale reliably without requiring sequential recurrence. Table 16: Sensitivity to extended forecasting horizons (Lp=192L_p=192). Bold: best; underline: second best. LfL_f Method ETTh1 ETTh2 Weather Electricity Traffic† CRPS (×10−2× 10^-2, ↓ ) 336 Raw 17.77 ± 0.23 13.07 ± 0.33 3.98 ± 0.05 4.43 ± 0.09 10.36 ± 0.20 Fourier 18.07 ± 0.49 12.91 ± 0.19 3.82 ± 0.02 4.51 ± 0.10 10.41 ± 0.15 DCT 17.81 ± 0.67 12.65 ± 0.31 3.82 ± 0.02 4.49 ± 0.09 10.45 ± 0.15 PCA 17.93 ± 0.49 12.51 ± 0.30 3.83 ± 0.03 4.49 ± 0.08 10.29 ± 0.16 MercerFlow 17.82 ± 0.49 12.66 ± 0.36 3.82 ± 0.03 4.52 ± 0.10 10.41 ± 0.21 720 Raw 18.61 ± 0.43 11.28 ± 0.27 5.13 ± 0.08 4.98 ± 0.36 12.06 ± 0.04 Fourier 18.48 ± 0.58 10.35 ± 0.13 4.78 ± 0.07 4.35 ± 0.12 12.03 ± 0.07 DCT 18.52 ± 0.77 11.52 ± 0.55 4.79 ± 0.03 4.33 ± 0.13 11.94 ± 0.07 PCA 18.77 ± 0.87 10.99 ± 0.54 4.76 ± 0.05 4.43 ± 0.11 11.99 ± 0.13 MercerFlow 18.86 ± 0.76 11.22 ± 0.54 4.79 ± 0.07 4.58 ± 0.21 12.06 ± 0.10 MAE (↓ ) 336 Raw 2.18 ± 0.03 5.19 ± 0.12 13.25 ± 0.29 187.6 ± 3.9 4.70 ± 0.10 Fourier 2.20 ± 0.05 5.12 ± 0.06 13.11 ± 0.14 190.5 ± 4.0 4.70 ± 0.10 DCT 2.18 ± 0.09 5.00 ± 0.12 13.06 ± 0.18 190.6 ± 3.7 4.70 ± 0.10 PCA 2.20 ± 0.06 4.94 ± 0.13 13.08 ± 0.19 190.5 ± 3.7 4.60 ± 0.10 MercerFlow 2.19 ± 0.06 5.01 ± 0.15 13.07 ± 0.15 191.3 ± 4.2 4.70 ± 0.10 720 Raw 2.38 ± 0.07 4.50 ± 0.09 15.37 ± 0.38 206.9 ± 15.6 5.30 ± 0.00 Fourier 2.40 ± 0.10 4.20 ± 0.07 14.33 ± 0.21 176.2 ± 6.6 5.30 ± 0.00 DCT 2.42 ± 0.13 4.65 ± 0.23 14.60 ± 0.25 175.7 ± 5.6 5.20 ± 0.00 PCA 2.46 ± 0.14 4.47 ± 0.21 14.51 ± 0.29 182.0 ± 5.3 5.30 ± 0.01 MercerFlow 2.48 ± 0.13 4.54 ± 0.24 14.63 ± 0.26 189.1 ± 9.6 5.30 ± 0.01 † Traffic MAE is reported in ×10−3× 10^-3. Appendix E Optimizer Dynamics and Stability Increasing the AdamW learning rate from the default η=10−3η=10^-3 to η=10−1η=10^-1 subjects the first-layer curvature to extreme stress under the shared OU prior. Table 17 reports test CRPS at a matched mid-run evaluation checkpoint across all five benchmarks (seeds 4242–4444). To ensure fair comparison, the shared evaluation checkpoint corresponds to the median of validation-best epochs among decorrelating codecs for each dataset (ETTh1 at epoch 4040, ETTh2 at 4040, Weather at 1010, Electricity at 2020, and Traffic at 3030). Table 17: Test CRPS (×10−2× 10^-2, ↓ ) at AdamW learning rate η=10−1η=10^-1 (default is 10−310^-3), mean ± std over seeds 4242–4444. Evaluated at matched mid-run EMA checkpoints E∗E^*. Bold: within one std of the best. Dataset Mercer (Ours) Fourier DCT PCA Raw ETTh1 17.99±0.4417.99± 0.44 18.17±0.5518.17± 0.55 18.51±0.3318.51± 0.33 17.81±0.0017.81± 0.00 122.42±80.10122.42± 80.10 ETTh2 13.05±0.5613.05± 0.56 13.27±0.0613.27± 0.06 14.00±1.1314.00± 1.13 12.92±0.2012.92± 0.20 15.18±0.6915.18± 0.69 Weather 1.70±0.011.70± 0.01 1.71±0.011.71± 0.01 1.71±0.011.71± 0.01 1.70±0.011.70± 0.01 1.71±0.061.71± 0.06 Electricity 4.34±0.064.34± 0.06 4.33±0.064.33± 0.06 4.34±0.044.34± 0.04 4.39±0.054.39± 0.05 6.53±0.356.53± 0.35 Traffic 10.65±0.3610.65± 0.36 10.80±0.8610.80± 0.86 11.16±0.2611.16± 0.26 10.65±0.4110.65± 0.41 24.62±7.1724.62± 7.17 Under this severe learning rate regime, the unrotated Raw baseline—whose first-layer condition number κ(DAD)κ(DAD) is two to four orders of magnitude larger than that of the orthonormal codecs—catastrophically collapses on the harder benchmarks (CRPS explodes to 1.221.22 on ETTh1, 0.0650.065 on Electricity, and 0.2460.246 on Traffic). In stark contrast, all decorrelating codecs (MercerFlow, Fourier, DCT, PCA) preserve stable descent trajectories, avoiding numerical instability and remaining close to their default-learning-rate performance. To visualize optimization stability, Figure 3 plots the training loss trajectories at an elevated learning rate (η=10−2η=10^-2). MercerFlow converges rapidly and strictly monotonically. In contrast, the unrotated baseline (Raw, labelled RawFlow in the legend) exhibits pronounced optimisation friction and overt instability (e.g., a severe gradient-induced loss spike around epoch 6060 on ETTh1), visually confirming the critical role of first-layer preconditioning. Figure 3: Training loss trajectories on ETTh1 and ETTh2 at an elevated learning rate (η=10−2η=10^-2). MercerFlow (solid blue) converges faster and strictly monotonically. In contrast, the unrotated baseline (RawFlow, dashed orange) suffers from optimization friction and visible instability (e.g., the loss spike around epoch 6060 on ETTh1), demonstrating the empirical benefits of first-layer preconditioning. Appendix F Extended Ablations F.1 Synthetic Ornstein–Uhlenbeck Comparison Evaluating models on synthetic Ornstein–Uhlenbeck (OU) sample paths provides a strictly controlled testbed with a known data-generating kernel. Under this setting, matched and mismatched linear codecs are analytically well-defined, allowing us to evaluate the impact of first-layer conditioning in isolation from real-world confounding factors. Synthetic trajectories are generated with context length Lp=32L_p=32, forecast horizon Lf=32L_f=32 (total window L=64L=64), and stationary length scale ℓ=32 =32 (corresponding to one e-fold decay across the boundary between context and horizon). The observation-space GP prior and the Mercer basis share this identical ground-truth kernel. Fourier, DCT, and the unrotated Raw baseline employ the exact same AdaLN MLP backbone (width 6464, 22 hidden layers) and the same conditional GP prior; only the linear coordinate transformation Φ∈ℝL×L ^L× L differs. PCA is omitted in this setting as the ground-truth kernel is known analytically. All models are trained with AdamW at η=10−3η=10^-3 for 2,0002,000 steps (batch size 6464, random seeds 00–22). Table 18 reports the Jacobi condition number κ(DAD)κ(DAD) evaluated on the centred prior population (Φ⊤KΦ K ) and on empirical conditional CFM interpolants ctc_t, along with test CRPS (2020 Euler integration steps, 3232 trajectories, 256256 evaluation windows). Matched OU Prior. On the centred prior covariance, MercerFlow analytically achieves optimal conditioning (κ(DAD)=1.00κ(DAD)=1.00). DCT closely follows (κ(DAD)=1.82κ(DAD)=1.82), Fourier sits in the tens (24.924.9), and the unrotated Raw baseline suffers from severe ill-conditioning (κ(DAD)=1,501.0κ(DAD)=1,501.0). On empirical conditional interpolants, conditioning degrades across all codecs due to non-zero context conditioning, but preserves the identical relative ordering: Mercer (98.598.5) and DCT (105.0105.0) remain around 10210^2, Fourier sits at 654.0654.0, and Raw reaches 3,515.03,515.0. Under AdamW, all four codecs tie on test CRPS within 0.2%0.2\% (DCT 0.6062±0.00240.6062± 0.0024, Mercer 0.6063±0.00250.6063± 0.0025, Fourier 0.6064±0.00230.6064± 0.0023, Raw 0.6069±0.00300.6069± 0.0030). This flat scoreboard reflects the synthetic “easy” regime: because the prior matches the data-generating process exactly, the prior conditional mean already explains almost all predictive variance, leaving minimal residual dynamics for the flow to resolve. Nevertheless, the geometric separation remains stark in the condition number columns. Mismatched Prior Length Scale (Wrong-ℓ Stress Test). To examine sensitivity to kernel misspecification, we train a mismatched MercerFlow arm where the GP model and Mercer basis assume ℓ=320 =320 (10×10× larger than the ground-truth data length scale ℓ=32 =32). This mismatch increases the prior κ(DAD)κ(DAD) from 1.001.00 to 1.711.71 and the interpolant κ(DAD)κ(DAD) from 98.598.5 to 192.0192.0, resulting in a mild CRPS degradation (0.6063→0.6125±0.00570.6063→ 0.6125± 0.0057). Crucially, this performance gap (0.00620.0062) remains within seed standard deviation, confirming that MercerFlow maintains stable optimization even under moderately misspecified prior parameters. SGD Equivariance Control. Training with Nesterov SGD (momentum 0.90.9, learning rate 10−310^-3, matched step budget) completely eliminates all performance discrepancies: Mercer, Fourier, and Raw achieve identical test CRPS (0.6033±0.00340.6033± 0.0034), with a relative gap from Raw of exactly 0.00%0.00\%. This empirical result provides definitive proof of Proposition 1, demonstrating that coordinate representations have zero effect under rotation-equivariant optimization. Table 18: Synthetic OU comparison (context/horizon 32/3232/32, data length scale ℓ=32 =32). Matched codecs share the OU prior at ℓ=32 =32; the wrong-ℓ arm uses a Mercer/GP model at ℓ=320 =320. AdaLN MLP (6464 hidden, 22 layers), AdamW η=10−3η=10^-3, 2,0002,000 steps, seeds 00–22 (mean ± std). κ(DAD)κ(DAD) prior is computed on the population covariance Φ⊤KΦ K ; interpolant κ(DAD)κ(DAD) is empirical on conditional CFM states ctc_t. SGD control: Nesterov momentum 0.90.9, matched compute budget; gaps are 100×(CRPS/CRPSRaw−1)100×(CRPS/CRPS_Raw-1) (mean over seeds). Bold indicates values within one standard deviation of the best. Codec / arm κ(DAD)κ(DAD) prior κ(DAD)κ(DAD) interp. Test CRPS (↓ ) SGD gap (%) Mercer (ℓ=32 =32) 1.001.00 98.598.5 0.6063±0.00250.6063± 0.0025 0.000.00 Fourier (ℓ=32 =32) 24.9024.90 654.0654.0 0.6064±0.00230.6064± 0.0023 0.000.00 DCT (ℓ=32 =32) 1.821.82 105.0105.0 0.6062±0.00240.6062± 0.0024 0.000.00 Raw (ℓ=32 =32) 1,501.01,501.0 3,515.03,515.0 0.6069±0.00300.6069± 0.0030 0.000.00 Mercer (ℓ=320 =320) 1.711.71 192.0192.0 0.6125±0.00570.6125± 0.0057 — F.2 Sampling Efficiency: Flow Matching vs. Diffusion Conditional flow matching generates straight probability paths by pairing source prior draws directly with data targets, avoiding the stochastic Brownian simulation inherent to diffusion. Table 19 benchmarks MercerFlow against an equivalent score-based diffusion model formulated in the identical Mercer latent space (MercerDiff) across varying Numbers of Function Evaluations (NFE) at forecast horizon Lf=192L_f=192. While MercerDiff exhibits severe discretization errors at low step budgets and remains suboptimal even at NFE=128NFE=128 (0.470.47 on ETTh1), MercerFlow fully saturates its predictive distribution within merely 44–1010 Euler integration steps (0.1990.199), providing a massive practical acceleration during inference. Table 19: Inference sampling efficiency: Test CRPS (↓ , mean ± std) across varying Numbers of Function Evaluations (NFE) at horizon Lf=192L_f=192. Due to straight probability paths, MercerFlow saturates within 44–1010 steps, whereas diffusion (MercerDiff) suffers from severe discretization error. Bold: best per column. NFE (Euler Steps) Dataset Method 4 10 20 50 128 ETTh1 MercerDiff 19.03±0.0719.03± 0.07 1.04±0.061.04± 0.06 0.56±0.070.56± 0.07 0.51±0.060.51± 0.06 0.47±0.030.47± 0.03 MercerFlow 0.200±0.0170.200± 0.017 0.199±0.0170.199± 0.017 0.199±0.0170.199± 0.017 — — ETTh2 MercerDiff 7.78±0.027.78± 0.02 0.47±0.010.47± 0.01 0.211±0.0050.211± 0.005 0.176±0.0050.176± 0.005 0.177±0.0030.177± 0.003 MercerFlow 0.137±0.0050.137± 0.005 0.136±0.0060.136± 0.006 0.136±0.0070.136± 0.007 — — F.3 Sensitivity to Prior Diagonal Jitter The regularising diagonal jitter εIL I_L added to the Gaussian process covariance matrix K serves two distinct purposes: it ensures numerical stability during matrix inversion in the conditional prior (Equation 3) and sets the diagonal regularisation floor for the Mercer eigenspace. A minimal jitter level (ε≤10−3 ≤ 10^-3) preserves exact alignment with the continuous analytical kernel, whereas excessive jitter forces the prior covariance toward isotropic white noise (K≈ILK≈ I_L). Table 20 systematically evaluates the sensitivity of first-layer conditioning κ(DAD)κ(DAD) (measured on conditional interpolants tc_t) and test CRPS across jitter values ε∈[10−4,1.0] ∈[10^-4,1.0] on ETTh1 and ETTh2. Stability Under Standard Regularisation (ε≤10−2 ≤ 10^-2). For standard jitter scales (ε∈[10−4,10−2] ∈[10^-4,10^-2]), MercerFlow maintains tightly bounded, well-conditioned first-layer geometry (κ(DAD)≤30.6κ(DAD)≤ 30.6 on ETTh1 and ≤95.3≤ 95.3 on ETTh2), yielding stable predictive performance (CRPS 0.1860.186–0.1880.188 on ETTh1). Conversely, diagonal jitter fails to resolve the conditioning pathology of the unrotated Raw baseline, whose condition number remains several orders of magnitude higher (κ(DAD)∼105κ(DAD) 10^5–10610^6) across all evaluated jitter scales. On the simple ETTh2 benchmark, where the GP prior mean accounts for the dominant dynamics, CRPS remains virtually invariant to the choice of codec (0.1310.131–0.1350.135). Breakdown Under Extreme Jitter (ε≥0.1 ≥ 0.1). At extreme jitter levels (ε≥0.1 ≥ 0.1), the kernel’s continuous temporal correlation is drowned out by artificial diagonal dampening (K→ILK→ I_L). This distorts the kernel’s eigenspace, driving Mercer’s condition number up to 3.64×1033.64× 10^3 on ETTh2. More critically, because the prior degenerates toward uninformative white noise, the flow loses its data-matched starting state, causing pronounced CRPS degradation on the harder ETTh1 benchmark (0.188→0.2550.188→ 0.255 for Mercer and 0.192→0.2560.192→ 0.256 for Raw). These results show that MercerFlow’s effectiveness stems from the kernel’s authentic analytical structure rather than artificial diagonal conditioning. Table 20: Sensitivity to GP prior diagonal jitter ε on ETTh1 and ETTh2 (L=384L=384). κ(DAD)κ(DAD) is evaluated on conditional interpolants tc_t. Standard regularisation (ε≤10−2 ≤ 10^-2) preserves well-conditioned representations; extreme jitter (ε≥0.1 ≥ 0.1) collapses the prior toward white noise (K≈ILK≈ I_L), degrading predictive accuracy. Bold: best per column within each dataset block. Diagonal Jitter ε Metric Dataset Codec 10−410^-4 10−310^-3 10−210^-2 10−110^-1 1.01.0 κ(DAD)κ(DAD) (tc_t, ↓ ) ETTh1 Mercer 23.55 30.64 29.37 47.03 334.67 Raw 2.01×1052.01× 10^5 3.02×1053.02× 10^5 1.85×1051.85× 10^5 1.66×1051.66× 10^5 1.44×1051.44× 10^5 ETTh2 Mercer 62.14 74.75 95.27 524.26 3.64×3.64× 10^3 Raw 4.07×1054.07× 10^5 1.22×1061.22× 10^6 1.20×1061.20× 10^6 1.17×1061.17× 10^6 1.01×1061.01× 10^6 CRPS (↓ ) ETTh1 Mercer 0.188 0.186 0.187 0.195 0.255 Raw 0.192 0.187 0.187 0.196 0.256 ETTh2 Mercer 0.133 0.132 0.135 0.132 0.136 Raw 0.131 0.135 0.134 0.131 0.136