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.
We propose a closed-form spectral framework for relative log-density estimation in linearly parameterized probabilistic models, including unnormalized and conditional models. This is achieved by representing the Kullback-Leibler (KL) divergence as an integral of weighted chi-squared divergences, converting KL estimation into a family of least-squares problems. We derive an explicit spectral formula based only on first- and second-order feature moments, yielding closed-form estimators of both divergences and log-density potentials for fixed features. The framework extends to a broad class of f-divergences and can be combined with kernelization or feature learning with neural networks. We prove convergence guarantees for the resulting estimators and empirically compare them on synthetic data with optimization-based variational formulations, including logistic and softmax regression for normalized conditional models.
Francis Bach
SIERRA · Inria - Ecole Normale Supérieure · PSL Research University
We study ridge-regularized log-density-ratio estimation in the Gaussian location model with a common covariance matrix. By affine invariance, the model is written as q ∼ N(0, I), p ∼ N(Δ, I), with linear features, where Δ is a mean vector. The variational estimator is the empirical Kullback-Leibler (KL) log-normalized fit with a squared L2-penalty on its nonconstant coefficient, and the spectral estimator recently introduced in [1] replaces a single variational problem by a continuum of ridge-regularized least-squares problems. We derive high-dimensional deterministic asymptotic equivalents when the numbers of observations and dimension tend to infinity with fixed ratios. The regularized variational limit is characterized by a scalar entropy minimization problem derived from the convex-Gaussian-min-max theorem (CGMT), while the regularized spectral limit follows from deterministic equivalents for resolvents of weighted sums of two independent Gaussian sample covariance matrices. We use these formulas to compare population risks, with experiments focused on fixed-signal aspect-ratio sweeps and optimized regularization. Our conclusion is that with many observations, under the criteria and asymptotic regimes analyzed here, the well-specified variational estimator has the smaller risk, while with fewer observations, the spectral estimator is favored because its covariance-based construction has lower variance. We also study how a nuclear penalty can be used and partially analyzed to perform feature learning.
Francis Bach
SIERRA · Inria - Ecole Normale Sup´erieure PSL Research University
Density estimation is a fundamental problem in statistics and machine learning. In this work, we introduce BayesNDE, a neural density estimator based on Bayesian generative modeling. BayesNDE learns a Bayesian generative model and evaluates its density without requiring invertible networks or Jacobian-determinant computation. For each observation, it infers a sample-specific latent posterior to construct an adaptive proposal that focuses computation on regions contributing most to its density. Bridge sampling then combines samples from this proposal with separate posterior samples to estimate the density. Experiments on nonlinear and multimodal synthetic datasets show improved estimation of density values and better recovery of the density structure compared to the state-of-the-art neural density estimators. Applications to real-world datasets further demonstrate improved anomaly detection. Together, these results highlight BayesNDE as a flexible and effective neural density estimator, demonstrating how posterior inference can turn generative models into tools for density estimation. The code and tutorials are available at https://github.com/liuq-lab/BayesNDE.
Chenglin Li, Qiao Liu
Department of Biostatistics, Yale University New Haven, Connecticut, USA