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.
College of Computer and Information Technology, China Three Gorges University, China · City University of New York, USA · University of California, Berkeley, USA +1
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