Conditional density estimation targets the full distribution of a response given covariates, as required, for example, for per-galaxy photometric redshifts. We develop a scalable Bayesian estimator based on the logistic Gaussian process. The log conditional density has a separable covariance: a Matérn kernel along the response, represented in a truncated Fourier basis on a circle, and a covariate kernel represented by Nyström features, which accommodate non-stationary kernels with input-dependent amplitudes and length scales. Instead of a Laplace or variational approximation, we sample the latent field of this finite-feature model. Given the hyperparameters, its posterior is strongly log-concave with a uniformly bounded Hessian, and we draw from it by simulating kinetic Langevin dynamics with symmetric minibatch splitting in Kronecker-whitened coordinates. Marginal-likelihood gradients follow from Fisher's identity as posterior expectations. Under the conditions of our analysis their bias is controlled by the sampler's step size and run length, and the predictive averages over the non-Gaussian latent posterior instead of a Gaussian around its mode. On photometric-redshift benchmarks with up to 3.9 million training observations, trained on a single GPU, the estimator is competitive with state-of-the-art tabular foundation models on density and calibration metrics.
Figures & tables
model
NLL
RMSE
CRPS
CDE
PIT-KS
N=400
LGP (ours)
−0.0256
0.7262
0.3635
−1.179
0.0471
s.d.
0.0388
0.0072
0.0074
0.049
0.0268
TabPFN-3.5
−0.0091
0.7085
0.3496
−1.129
0.0428
s.d.
0.0317
0.0021
0.0023
0.047
0.0258
TabICL-2
+0.2101
0.7079
0.3511
+0.055
0.0448
Table 1: Bimodal sinusoid: test metrics on 2,000 test rows of the vanilla LGP (RBF kernels, Nyström features, three learned hyperparameters), TabPFN-3.5 and TabICL-2 with the training rows as context, Flow-Spline tuned by the protocol of Izbicki and Rodrigues (2026) , the CRPS ensemble (batch 32 , early stopping on 10% of the training rows with patience 200 epochs), and the true density. All rows are averaged over four seeds, each a different split of one sample into training and test rows (standard deviations in the rows below). A PIT-KS of about 0.02 is the sampling noise of 2,000 rows. Bold: best of the five models.
benchmark (train / test)
model
NLL
RMSE
CRPS
CDE
PIT-KS
NMAD
outliers
Happy A → B
LGP (ours)
−2.1557
0.0506
0.0203
−12.86
0.0054
0.0186
0.86%
( 74,950 / 74,900 )
s.d.
0.0009
0.00008
0.00002
0.02
0.0003
0.00003
0.01%
TabPFN-3.5
−2.1632
0.0493
0.0200
−12.90
0.0040
0.0185
0.81%
TabICL-2
−1.9958
0.0497
0.0201
−8.79
0.0068
0.0185
0.83%
Flow-Spline
−2.1093
0.0522
0.0208
−12.04
0.0220
0.0191
0.91%
s.d.
0.0049
0.00015
0.00003
0.16
0.0101
0.00007
0.01%
Table 2: Test metrics on the four photometric-redshift benchmarks (NLL and CDE loss: lower is better; RMSE and CRPS in redshift units; PIT-KS: Kolmogorov–Smirnov distance of the PIT from uniform). Bold: best of the five models (all tied values at the displayed precision). Rows followed by an s.d. row are means over four seeds (five runs for ‡ ), with the standard deviation (s.d.) in that row. The seeds change the validation rows and each model’s initialisation (and Flow-Spline’s search), and on SDSS and DESI also the train/test split; TabPFN-3.5 and TabICL-2 on Happy and Teddy do not depend on the seed. TabICL-2 and Flow-Spline on DESI are single runs on the seed- 0 split, on which the LGP’s NLL is −2.3795 . ‡ Mean over five runs with identical settings on that split, two of which diverge (Appendix I ). ∗ A context of 500,000 training rows. † Flow-Spline’s exact density on DESI has narrow spikes that dominate ∫p2 ; on the benchmark’s 200 -point grid its CDE loss is −13.95 .
network
Jw
sweeps
batches
time
Happy , Teddy
3×32
6
2+4
10
6 min
SDSS
6×64
6
1+2
20
16 min
DESI
10×128
10
1+2
40
2.2–2.5 h
Table 3: GP configuration and cost on one RTX 5090 (32 GB), scoring included.
Appendix figures & tables6 assets
Supplementary material from the paper’s appendix.
Appendix
own predictive
sampled pred.
σ2
(ℓx,ℓy)
sampling
−0.192(0.014)
−0.192(0.016)
59 – 78
(0.13 – 0.15,0.11 – 0.12)
MAP
+1.577(0.225)
+1.308(0.166)
3983 – 5062
(0.027,0.010)
Laplace
−0.076(0.013)
−0.192(0.016)
75 – 94
(0.14 – 0.15,0.11 – 0.12)
Appendix
Table 4: Learning ϑ on the sinusoid ( N=2,000 , J=100 , K=50 ): mean (standard deviation) of the test NLL over four seeds, by each method’s own predictive and by the sampled predictive at its learned ϑ , and the learned variance and length scales. The seed standard deviations are neither a convergence certificate nor an estimate of the predictive’s Monte Carlo error.
Figure 1: The vanilla logistic GP on the bimodal sinusoid with N=400 training samples (orange points). The floor heatmap shows the posterior predictive density; purple profiles show the predictive at x=−0.7,0,0.7 . Gray basis curves on the x=−1 boundary illustrate cos(πkT(y)) and −sin(πkT(y)) for k=1,2,4,8,16 , evaluated using the fixed training-only KDE-CDF transform T . Each sine/cosine pair shares a display baseline. The basis curves are scaled and vertically offset for readability, and their heights are not on the density axis.
NLL per station
Δ NLL vs SLGP
model
NLL
CRPS
RMSE
Bern
Grimsel
Meiringen
paper ℓ
CV ℓ
SLGP, paper’s ℓ , MAP
3.4000
4.251
7.39
3.371
3.361
3.469
—
+0.016±0.009
SLGP, paper’s ℓ , NUTS
3.4007
4.251
7.39
3.371
3.361
3.470
+0.001±0.001
+0.017±0.008
SLGP, CV ℓ , MAP
3.3838
4.216
7.34
3.366
3.355
3.431
−0.016±0.009
—
SLGP, CV ℓ , NUTS
3.3841
4.217
7.34
3.366
3.356
3.431
−0.016±0.009
+0.000±0.000
LGP, stationary (5 seeds)
3.3812 ( 0.0001 )
4.200
7.32
3.360
3.347
3.436
−0.019±0.014
−0.003±0.011
Appendix
Table 5: Swiss temperatures, the 3 Bern test stations ( 1,095 days, ∘ C). NLL and CRPS: mean over the test days; NLL per station; Δ : paired NLL difference to SLGP with the paper’s length scales (MAP) and to SLGP with station-CV length scales (MAP), ± the standard error of a block bootstrap over the 52 weeks of the year (the days of a week at the three stations form a block, since daily temperatures are autocorrelated). Rows over several runs: mean over the runs, standard deviation in parentheses. SLGP: 500 basis functions, Matérn- 5/2 ; paper’s ℓ : the saved fits of Gautier and Ginsbourger (2026) (its package build); CV ℓ : fitted with CRAN SLGP 2.0.0; LGP: our model. Bold: the lowest value in each column except the Δ NLL columns.
Figure 2: Predictive densities at the three Bern test stations over a histogram of the 365 observed daily means of 2019, with each model’s NLL at the station. Each model predicts one density per station, from its coordinates.
model
NLL
excess
Hbin2
CRPS
PIT-KS
Δ NLL vs LGP
time
seeds
N=400
LGP, stationary
−0.020
0.233
0.0621
0.3608
0.045
—
6 min
4
SLGP, package defaults, MAP
+0.619
0.873
0.2903
0.3771
0.087
+0.640±0.081
8 s
4
SLGP, tuned, MAP
−0.037
0.217
0.0485
0.3658
0.039
−0.016±0.040
40 s
4
SLGP, tuned, Laplace
+0.632
0.886
0.3187
0.6367
0.144
+0.652±0.195
40 s
4
SLGP, tuned, NUTS ∗
−0.021
0.232
0.0524
0.3712
0.043
−0.008±0.041
6.9 h
3
Appendix
Table 6: Bimodal sinusoid, four seeds of 2,000 test rows each: mean over the seeds of the NLL, the excess NLL over the true density on the same rows, the binned squared Hellinger distance Hbin2 to it, CRPS and PIT-KS; Δ : paired NLL difference to our model, mean ± standard deviation over the seeds (positive: worse than ours); time per fit, SLGP without its search; seeds: the number of seeds averaged. Every SLGP fit uses the exact basis scheme; tuned: length scales and σ2 selected on held-out training rows, in a search space designed after test-scored probe fits on seed 42 (see Protocol), so the comparison is exploratory. ∗ Seeds 42 – 44 ; on seed 45 NUTS did not converge (see text).
Figure 3: Happy A → B: test NLL on 10,000 galaxies of sample B against the number of training galaxies (left) and against the fit time (right; along each line n increases; our model on one A100 including compilation, SLGP on two CPU cores, excluding the 26 -minute search that selected its hyperparameters). “Tuned” SLGP uses the length scales and σ2 selected on 1,000 galaxies.