Neural network-based predictive modeling with high-dimensional structured Gaussian targets requires an efficient and numerically stable, yet expressive approximation of the covariance matrix. We propose SCORE: a scalable framework, combining scoring rule training with an expressive covariance approximation learned in spectral space. For d-dimensional data, the learning task is decomposed into learning the marginal distributions and learning a structured correlation matrix, which enables dense dependencies with linear storage and O(dlogd) cost. We utilize the closed form Gaussian kernel score for training, which remains defined even for degenerate covariances and admits bounded gradients during optimization. We characterize kernel scores under invertible transforms and prove exact invariance under unitary transforms. At population level, our two-level objective recovers the true marginals and projects the target correlation onto the representable class; finite-sample PAC bounds show that the errors of the two stages enter additively. We evaluate our model on a variety of tasks with a commonly assumed Gaussian domain: Time-series forecasting, monocular depth estimation, and spatial weather prediction, showing improved performance at lower computational cost.
Figures & tables
Figure 1: Schematic overview of SCORE . Stage 1 fits the mean μ and the marginal variances τ (orange) and standardizes the target, leaving only its correlation structure. Stage 2 models this correlation by a normalized spectrum sk⋆ in Fourier coordinates, whose back-transform F∗ is the correlation core R⋆ (purple). Rescaling by the frozen marginals yields the predictive covariance ΣP (green).
Params
Score eval.
Sampling
Σ^ij
Chol
d(d+1)/2
O(d2)
O(d2)
O(d)
Diag
d
O(d)
O(d)
O(1)
LorD
d+dr
O(dr2+r3)
O(dr)
O(r)
SCORE (Ours)
2d
O(dlogd)
O(dlogd)
O(d)
Table 1: Cost of covariance parameterizations: number of covariance parameters, one log-score evaluation, one sample, and estimation of one covariance entry Σ^ij .
Figure 2: Toy example. We draw a signal y∼N(0,Σpr) with circulant Σpr and impose a cut-off on the eigenvalues, i.e., λi>0 iff ∣i∣≤B . Now, we observe a noisy subsample x=Sy+η with observational noise η∼N(0,σ2In) , where S∈Rn×d selects every (d/n) -th coordinate, and train a neural network to reconstruct y given x . The posterior has a closed form with rank(Σpost)=min(2B+1,d) . By varying the band limit B , we interpolate between a low-rank and a higher-rank target covariance matrix. The figure shows the estimated mean (colored dashed line) and standard deviation (shaded colored area), as well as the true posterior standard deviation (shaded grey area). In the low-rank case, Slog fails to recover the posterior, while Sk does. More details are given in Appendix F .
Figure 3: SCORE predictions. Left: Mean prediction and standard deviation of SCORE on the univariate ETTh2 dataset and a comparison with LorD regarding similarity to the empirical correlation structure. Right: observed temperature across stations and predictions and correlations using SCORE and LorD.
n
MSE
CRPS
ES
VS
167
150
140
156
139
Table 4: Number of wins when training the same predictive model with the Gaussian kernel score instead of the NLL. Bold : p<0.05 under the clustered sign-flip test, Holm-Bonferroni corrected over the four metrics.
Figure 4: Performance improvement (difference in nats per coordinate for NLL, percentage improvement otherwise) vs. increase in compute, relative to the Diag (median over seeds and data sets).
Figure 5: Different calibration measures for LorD, SB and SCORE , using the kernel score across all datasets. The dashed diagonal line indicates optimal calibration.
Appendix figures & tables11 assets
Supplementary material from the paper’s appendix.
Appendix
Figure 6: Realizations of the signal y for different band limits B∈{4,8,16,32} (upper) and corresponding eigenvalues of the posterior covariance Σpost (lower).
B=4
B=8
B=16
B=32
Sk
0.00357±3.7⋅10−6
0.108±0.014
0.618±0.011
0.584±0.010
Slog ( ε=10−6 )
1.15±0.26
1.43±0.11
2.07±0.22
1.95±0.10
Slog
0.953±0.13
1.33±0.12
2.08±0.22
1.99±0.26
Appendix
Table 5: Frobenius distance ∥Σϕ−Σpost∥F between the learned covariance Σϕ and the true posterior covariance Σpost across different band limits B . The results are averaged across five seeds with the best model in bold.
Figure 7: True posterior and different estimates for a selected realization and varying band limit B∈{4,8,16,32} . The ranks of the posterior covariance are identical to those in Figure 6 .
Figure 8: Gradient norm of the different loss functions with respect to the estimated covariance Σϕ , compared across different band limits B .
Data set
Variables
Resolution
Observations
Split
L
T
ETTh1, ETTh2
7
hourly
17,420
12/4/4 m
336
{96,192}
ETTm1, ETTm2
7
15 min
69,680
12/4/4 m
336
{96,192}
Weather
21
10 min
52,696
70/10/20
336
{96,192}
Electricity
321
hourly
26,304
70/10/20
336
{96,192}
Traffic
862
hourly
17,544
70/10/20
336
{96,192}
ILI
7
weekly
966
70/10/20
104
{24,36}
Appendix
Table 7: Time-series data sets. Splits are chronological. ETT uses 12/4/4 months, the other data sets use 70/10/20 %.
Figure 9: Comparison of relative improvement and computational cost against the diagonal Gaussian baseline for the time series datasets. All methods are trained on the Gaussian kernel score. Each point represents the median performance over five seeds for one task, i.e., uni-/multivariate and prediction horizon T .
Figure 10: Comparison of relative improvement and computational cost against the diagonal Gaussian baseline for the remaining datasets. All methods are trained on the Gaussian kernel score. Each point represents the median performance over five seeds.
Figure 11: Empirical CDFs of the PIT values of the marginal, location and scale pre-rank functions for the different predictive methods. A calibrated prediction lies on the diagonal. Results are aggregated over five seeds.
Figure 12: Empirical CDFs of the PIT values of the variogram pre-rank function with respect to different spatial (or temporal) lags h for the different predictive methods. A calibrated prediction lies on the diagonal. Results are aggregated over five seeds.
Figure 13: Relative performance of the low-rank method against ours with increasing rank r .
Figure 14: Performance of staged training (3:1 split) vs. first-stage only using the same computational budget.
Direct forecasting has become a standard paradigm for multivariate time-series forecasting because it predicts the full future horizon in a single pass. However, its training objective is often still decomposed into pointwise errors such as MSE. Such objectives provide stable supervision, but they do not explicitly preserve the structure of the future trajectory: temporal coherence within each variable and relational consistency across variables can both be weakened. We propose CoRe, a model-agnostic learning objective for direct multivariate forecasting. CoRe replaces pointwise supervision with two output-space constraints: a frequency coherence loss that aligns predicted and target spectra, and a low-rank relational graph loss that matches sampled pairwise differences in a target-derived PCA subspace. The resulting objective introduces no trainable parameters and can be applied to existing forecasting backbones by changing only the loss. Experiments on standard benchmarks show that CoRe improves strong baselines, compares favorably with recent forecasting objectives, and remains effective across different backbones, datasets, and hyperparameter settings overall consistently.
Xiaoyu Lin, Huiran Duan, Yining Liu +3
College of Computer and Information Technology, China Three Gorges University, China · City University of New York, USA · University of California, Berkeley, USA +1
This feature article provides an overview of the theoretical foundations for coVariance neural networks (VNNs), i.e., graph neural networks (GNNs) operating on covariance matrices as graphs. Covariance matrices are ubiquitous across domains, and hence, the deployment of GNNs often leverages graphs of pairwise statistical dependencies. Existing theoretical contributions on GNNs consider abstract graph representations and cannot accommodate the data-driven nuances associated with covariance matrices. This tutorial brings into focus various novel theoretical insights via mathematical analyses of VNNs that have broad signal processing implications, including: (i) a conceptual equivalence between VNNs and principal component analysis (PCA)-based information processing; (ii) refined stability bounds on predictive outcomes in the presence of finite sample-induced covariance matrix perturbations; and (iii) refined characterization of transferability of VNNs across multiscale datasets. The theoretical insights discussed herein provide the underlying principles and justification towards adopting VNNs over workhorse PCA-based learning pipelines, in applications where covariance matrices are useful descriptors of data structure. We also convey how impact of these foundational advances permeates to \textit{principled} designs and applications of learning methods across broad domains where covariance matrices emerge. Notably, we elucidate the conceptual insights facilitated by VNNs to the specific task of characterizing brain age gap for neurodegenerative conditions using neuroimaging datasets, a timely problem in computational neuroscience. Broader impacts to other application domains are discussed as well.
Saurabh Sihag, Andrea Cavallo, Elvin Isufi +2
Department of Electrical and Computer Engineering at the University at Albany, SUNY, Albany, NY · Delft University of Technology, Delft, The Netherlands · Department of Electrical and Computer Engineering at the University of Rochester, Rochester, NY +1
We present a theoretically grounded Gaussian process framework that leverages neural feature maps to construct expressive kernels. We show that the learned feature map can be interpreted as an optimal low-rank approximation to a Gram matrix derived from an implied RKHS, from which we establish consistency of the GP posterior. We further analyse the spectral properties of the induced kernels and introduce product feature-map kernels to address oversmoothing. This simple yet powerful approach enables fast, scalable, and accurate exact GP inference with minimal upfront work. The flexibility of kernel design supports seamless application to both regression and classification tasks across diverse data modalities, including tabular inputs and structured domains such as images. On benchmark datasets, this approach surpasses pre-existing methods in terms of accuracy and training and prediction efficiency.