Pith. sign in

REVIEW 3 major objections 4 minor 3 references

Distributional encoding for Gaussian process regression with qualitative inputs

T0 review · 3 major / 4 minor · reviewed 2026-08-07 · deepseek-v4-flash

Pith's one-line read Distributional encoding puts categorical inputs into Gaussian process regression by representing each category as the empirical distribution of the response at that level, using positive semidefinite kernels on distributions inside the…

desk verdict Useful extension of target encoding for GPs with categorical inputs, but target leakage in the encoding and LOO selection could bias the reported gains. read the letter →

arxiv 2506.04813 v1 pith:M6W3G6ZO submitted 2025-06-05 stat.ML cs.LG

classification stat.MLcs.LG MSC 62G0862M2068T05
keywords GaussianprocessregressioncategoricalinputsdistributionalencodingtargetWassersteindistancemaximummeandiscrepancykernelmethodsmixed-variableBayesianoptimization
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

Distributional encoding is proposed as a way to put categorical inputs into Gaussian process regression without learning a latent coordinate for each category. Each level of a categorical variable is replaced by the empirical distribution of the training responses observed at that level, and a positive semidefinite kernel on distributions, built from the 2-Wasserstein distance or maximum mean discrepancy, is multiplied into the usual product kernel. The paper argues this keeps all GP inference unchanged while matching or beating a latent-variable GP on the tested benchmarks and doing so at much lower computational cost. The same recipe is extended to classification, multi-task outputs, and auxiliary or low-fidelity data, where the encoding can be formed from a separate data source and even handles new levels without retraining.

What carries the argument

The carrying object is a positive semidefinite kernel on probability distributions, used in place of the categorical similarity matrix $[T_j]_{u_j,u'_j}$ in the tensorized GP kernel. For scalar outputs the paper exploits the result that the distance substitution kernel $k_{W_2}(P,Q)=\exp(-\gamma W_2^\beta(P,Q))$ is positive semidefinite for $\beta\in[0,2]$; for MMD, it uses $k_{MMD}(P,Q)=\exp(-\gamma \operatorname{MMD}^2(P,Q))$, which is positive semidefinite in any dimension and reduces to the energy distance for a particular base kernel. These kernels convert empirical distributions of the response into distances between categorical levels, so the standard GP machinery, including the Gram matrix, marginal likelihood, and prediction, runs unchanged. For multivariate outputs, the univariate obstruction is bypassed with the sliced Wasserstein distance, and for classification the paper uses histogram kernels based on $\chi^2$, total-variation, or Hellinger distances.

What would settle it

Rerun the material-design and engineering benchmarks with leave-one-out validation where, for each fold, the empirical distribution of the held-out level is rebuilt from the remaining training points only; if RRMSE degrades substantially relative to the paper's leave-one-out selection, the reported advantage is not robust to target leakage.

Watch

Extended reading notes

Core claim

On the paper's own terms, the central claim is that a categorical level can be represented by the full conditional distribution of the response, not just its mean, and that kernels on those distributions make this representation usable inside a product GP kernel. Specifically, for a level $l$ of input $u_t$, the encoding is $\hat P^Y_{t,l} = \frac{1}{N_{tl}}\sum_i \delta_{y^{(i)}} 1_{u^{(i)}_t = l}$, and the Gram matrix multiplies $\exp(-\gamma W_2^\beta(\hat P^Y_{t,u_t}, \hat P^Y_{t,u'_t}))$ or $\exp(-\gamma \operatorname{MMD}^2(\hat P^Y_{t,u_t}, \hat P^Y_{t,u'_t}))$ across categorical factors. The paper claims this weakly supervised encoding is on par with the latent-variable GP on a material-design benchmark, often better on engineering test functions in terms of relative root mean-squared error, while being much cheaper to fit, and that auxiliary data can be incorporated simply by enriching the empirical distributions.

Load-bearing premise

The load-bearing premise is that using the training response $y$ to build the input encoding does not leak target information into the kernel; if recomputing the per-level distribution without the held-out point changes the leave-one-out comparisons, the accuracy gains would be partly an artifact of leakage.

Editorial extensions

If this is right

  • Categorical inputs can be handled by standard GP toolkits without adding parameters per level, since only the kernel on the encoded distributions is introduced.
  • Auxiliary or low-fidelity response data can be folded into the encoding directly, so predictions for previously unseen levels are possible whenever auxiliary observations exist, without retraining.
  • Fast leave-one-out formulas can select among encodings, such as mean, mean and standard deviation, W2, and MMD, for each categorical input, making the choice of representation data-driven.
  • On the tested engineering and material benchmarks, the W2 and MMD encodings match the predictive accuracy of a latent-variable GP while reducing fitting time.
  • The same construction extends naturally to Bayesian optimization over mixed discrete and continuous spaces, because the encoded categorical inputs become compatible with continuous kernels.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • The paper does not analyze target leakage, but because the encoding is built from the same response the GP later predicts, recomputing each level's distribution without the held-out point during leave-one-out selection is a direct robustness check that could change the reported comparisons.
  • The approach likely gains most when per-level sample sizes are large enough to make empirical distributions informative; on small per-level samples, mean encoding may be more stable, so sample-size-aware shrinkage of the empirical distributions is a natural extension.
  • The framework points toward hybrid models: distributional encodings could initialize a latent-variable fit, combining the interpretability of direct distances with the flexibility of learned latent coordinates.
  • For categorical inputs with strong interactions with quantitative inputs, the paper's proposed sensitivity-analysis refinement is only sketched; a full empirical study on large datasets would test whether the pre-screening recovers the missing interaction terms.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 4 minor

Summary. The paper proposes 'distributional encoding' (DE) for Gaussian process regression with qualitative inputs: each level of a categorical variable is represented by the empirical distribution of the training response values observed at that level, and the GP kernel is built as a product kernel whose categorical factor is a positive semi-definite kernel on these empirical distributions (Wasserstein-2 or MMD/energy distance). The authors argue that this generalizes target encoding, is less costly than latent-variable GP (LVGP), and can be extended to classification, multi-task learning, and auxiliary/low-fidelity data. The empirical sections compare DE with LVGP, mean encoding, and mean/standard-deviation encoding on four engineering benchmark functions and a materials design problem, reporting RRMSE over 50 replications. The paper claims state-of-the-art predictive performance and stresses DE's lower computational cost relative to LVGP.

Significance. If the empirical claims hold, the paper offers a simple, computationally cheap alternative to LVGP for GP regression with categorical inputs, with useful extensions to multi-task and auxiliary-data settings. The manuscript is reproducible (R code supplied) and the experimental protocol is extensive, with 50 replications and several test cases. However, the central 'state-of-the-art' claim is not supported by the reported results, and the target-dependent construction of the kernel raises a leakage concern that is not addressed. The main methodological contribution—using distribution kernels on empirical conditional response distributions—is reasonable, but the current evidence is weakened by the leave-one-out selection issue and by the lack of a leakage-free validation protocol.

major comments (3)
  1. [Section 3 and Section 4.1 / Appendix B] The BestLOO protocol in Section 4.1 and Appendix B appears invalid because the fast leave-one-out formulas of Dubrule (1983) and Ginsbourger and Schärer (2025) assume a fixed Gram matrix, whereas the distributional encoding \hat P^Y_{t,l} is built from all training targets. When point i is removed, the empirical distribution of every other point in the same level changes, so many entries of K change. The paper does not state that encodings are recomputed without the held-out point, nor does it justify why the fixed-kernel LOO formulas remain applicable. The reported BestLOO RRMSE in Figure 6 can therefore be optimistically biased. Please recompute encodings within each LOO fold (or use nested cross-validation) and report whether the material-design results are stable under a leakage-free protocol.
  2. [Abstract and Section 4.1] The abstract's claim of 'state-of-the-art predictive performance' is contradicted by the paper's own results in Figure 5: LVGP has better RRMSE on Beam bending, mean encoding is best on Borehole, and W2/MMD encodings win on OTL and Piston. The material design experiment (Figure 6) shows DE on par with LVGP, not superior. The claim should be revised to something like 'competitive predictive performance at lower computational cost' and should explicitly acknowledge that no single method dominates across all test cases.
  3. [Section 2.3 and Section 3] Target encoding and distributional encoding are called 'weakly supervised' but they are constructed directly from the response y. This is a supervised construction, and because the GP kernel then depends on y, the posterior equations (1) are no longer a standard Bayesian update under a y-independent GP prior. The paper gives no theoretical justification or empirical safeguard (beyond the problematic LOO in the first comment) that this does not induce target leakage or overfitting. Please clarify what 'weakly supervised' means here, and provide at least an out-of-fold encoding comparison or a calibration/overfitting check to show that the reported test-set RRMSEs are not optimistic.
minor comments (4)
  1. [Definition 1] There is a typo: 'Wassertein' should be 'Wasserstein'. Also, in the estimator for \hat W_r, the symbol Q is used both for the number of quantile points and for one of the two distributions, which is confusing.
  2. [Section 3, Table 4] The illustrative table for distributional encoding appears broken: the column for the encoded input X3 is blank, so it does not actually show the distributional encoding values. Please fix the table formatting.
  3. [Section 2.3] The sentence 'In both cases, the resulting kernel will in fact be equivalent to the dirac kernel up to a multiplicative constant' is imprecise: when a Gaussian kernel is applied to one-hot vectors, the distance between different levels is a positive constant, leading to a kernel with equal off-diagonal entries, which is not exactly the Dirac kernel. Please specify the kernel family used for this equivalence.
  4. [Section 2.3] The text mentions 'a recent extension was proposed by Saves et al. (2023), with much better accuracy' but this method is not included in the benchmark comparisons. Either include it in the experiments or qualify the statement so it does not read as an unsupported comparison.

Circularity Check

0 steps flagged · score 0.0 of 10

No circularity: distributional encoding is a supervised feature map, but the GP predictions do not reduce to the encoding by construction; the fast-LOO selection raises a leakage risk that is a correctness concern, not a circular derivation.

full rationale

I walked the paper's derivation chain and found no step in which a claimed prediction or first-principles result reduces to its own inputs by construction. The encoding is defined from the training response, e.g. Eq. (4) of Section 3 with \hat P^Y_{t,l} = (1/N_{tl}) \sum_i \delta_{y_i} 1_{u_i=l}, and the Gram matrix then uses kernels between these empirical distributions. This is a supervised or weakly supervised feature map, but it is not a fitted parameter renamed as a prediction: the GP posterior mean is still obtained by standard conditioning on held-out targets, and no equation in the paper equates the predicted test value to the encoding itself. The positive-definiteness of the W2 and MMD kernels is imported from external cited theorems (Bachoc et al., 2017; Smola et al., 2007; Meunier et al., 2022), not from a self-citation, and the benchmarks compare against an external LVGP implementation on held-out test splits. The manuscript's limitation notes in Appendix A and the Conclusion concern small main effects and theoretical convergence, which are honest scope statements rather than circular moves. The fast leave-one-out selection in Section 4.1 and Appendix B is a genuine methodological risk: the fixed-kernel LOO formulas of Dubrule (1983) and Ginsbourger & Schärer (2025) are not exact when omitting a point changes the empirical distribution used in the kernel, so the 'BestLOO' criterion may be optimistic. But that is a statistical validity and implementation issue, not an equation-level circularity: the reported test RRMSE is still out-of-sample, and accepting the leakage concern would lower confidence in the empirical claim without making the derivation self-referential. There are no load-bearing self-citations and no imported uniqueness theorem. I therefore set the circularity score to 0.

Assumptions & free parameters 3 free parameters · 3 assumptions · 0 invented entities

The central claim depends on the PSD properties of known distribution kernels and on the domain assumption that the marginal response distribution per level captures the relevant behavior of a categorical input. The method also depends on the unstated choice of gamma, beta, and standard GP hyperparameters.

free parameters (3)
  • gamma (W2/MMD kernel lengthscale)
    Controls the scale of exp(-gamma*D^beta); the paper does not state whether it is fixed, estimated by maximum likelihood, or cross-validated, yet predictions depend on it.
  • beta (W2 kernel exponent)
    Allowed in [0,2] to maintain positive semi-definiteness, but no value or estimation method is given.
  • GP hyperparameters (lengthscales, nugget)
    Estimated by marginal likelihood via RobustGaSP. Standard for GP regression, but they are fitted to the same data and interact with the encoded inputs.
assumptions (3)
  • standard math PSD of the W2 distance kernel on R
    The method relies on Bachoc et al. (2017) that exp(-gamma W2^beta) is positive semi-definite on R for beta in [0,2].
  • standard math PSD of the MMD and energy distance kernel
    Relies on Song (2008) and Sejdinovic et al. (2013) that the Gaussian substitution of MMD squared with the energy distance kernel yields a PSD kernel in any dimension.
  • domain assumption Response distribution P(Y|u=l) characterizes a category's behavior
    Core encoding premise: using the marginal distribution of the target per level as the input representation is sufficient for GP similarity. Interactions with quantitative inputs are handled only in Appendix A as a workaround requiring large data.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Distributional encoding for Gaussian process regression with qualitative inputs." pith.science (2026). https://pith.science/paper/M6W3G6ZO

@misc{pith2026250604813,
  author       = {Pith},
  title        = {Pith review of: Distributional encoding for Gaussian process regression with qualitative inputs},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/M6W3G6ZO}},
  note         = {Machine review of arXiv:2506.04813}
}
read the original abstract

Gaussian Process (GP) regression is a popular and sample-efficient approach for many engineering applications, where observations are expensive to acquire, and is also a central ingredient of Bayesian optimization (BO), a highly prevailing method for the optimization of black-box functions. However, when all or some input variables are categorical, building a predictive and computationally efficient GP remains challenging. Starting from the naive target encoding idea, where the original categorical values are replaced with the mean of the target variable for that category, we propose a generalization based on distributional encoding (DE) which makes use of all samples of the target variable for a category. To handle this type of encoding inside the GP, we build upon recent results on characteristic kernels for probability distributions, based on the maximum mean discrepancy and the Wasserstein distance. We also discuss several extensions for classification, multi-task learning and incorporation or auxiliary information. Our approach is validated empirically, and we demonstrate state-of-the-art predictive performance on a variety of synthetic and real-world datasets. DE is naturally complementary to recent advances in BO over discrete and mixed-spaces.

Figures

Figures reproduced from arXiv: 2506.04813 by the authors.

Figure 1
Figure 1. Two-dimensional LVGP embedding versus mean/standard deviation encoding [PITH_FULL_IMAGE:figures/full_fig_p007_1.png] view at source ↗
Figure 2
Figure 2. Normalized distance matrix between levels for LVGP, mean and mean/standard [PITH_FULL_IMAGE:figures/full_fig_p007_2.png] view at source ↗
Figure 3
Figure 3. LVGP embedding versus W2 encoding for the beam bending test case [PITH_FULL_IMAGE:figures/full_fig_p010_3.png] view at source ↗
Figures from the paper (9 more)
Figure 4
Figure 4. Figure 4: Normalized distance matrix between levels for LVGP and [PITH_FULL_IMAGE:figures/full_fig_p010_4.png]
Figure 5
Figure 5. Figure 5: RRMSE on the four engineering test cases with 50 replications. [PITH_FULL_IMAGE:figures/full_fig_p016_5.png]
Figure 6
Figure 6. Figure 6: RRMSE on the material design test case with 50 replications. [PITH_FULL_IMAGE:figures/full_fig_p017_6.png]
Figure 7
Figure 7. Figure 7: RRMSE on the multi-output Borehole test case with 50 replications. Multi 1D 2D Output 1 0.215 (0.06) 0.21 (0.06) Output 2 0.221 (0.06) 0.215 (0.06) [PITH_FULL_IMAGE:figures/full_fig_p017_7.png]
Figure 8
Figure 8. Figure 8: RRMSE on the multi-fidelity Borehole test case with 50 replications. [PITH_FULL_IMAGE:figures/full_fig_p018_8.png]
Figure 9
Figure 9. Figure 9: f1 function: true function and training samples (left), output conditional prob￾ability for each level (right). Each level is represented by a different color. Before discussing a workaround, we can first try to detect if a qualitative input has a small main effect. In…
Figure 10
Figure 10. Figure 10: Computation time on the four engineering problems with 50 replications. [PITH_FULL_IMAGE:figures/full_fig_p025_10.png]
Figure 11
Figure 11. Figure 11: RRMSE versus computation time on the material design problem with 50 [PITH_FULL_IMAGE:figures/full_fig_p026_11.png]
Figure 12
Figure 12. Figure 12: RRMSE versus computation time on the material design problem with 50 [PITH_FULL_IMAGE:figures/full_fig_p026_12.png]

Discussion (0). Sign in to comment.

Reference graph

Works this paper leans on

3 extracted references · 3 canonical work pages

  1. [1]

    Compute first-order indices of all quantitative and qualitative inputs

  2. [2]

    > 5 − 10%), use the standard distributional encoding (b) Otherwise, compute its second-order interaction with all quantitative inputs

    For all qualitative inputs: (a) If it has a sufficiently large main effect (e.g. > 5 − 10%), use the standard distributional encoding (b) Otherwise, compute its second-order interaction with all quantitative inputs. For each significant interaction (e.g.> 5%), use distributional encoding on ˆP Y l,s = 1 Nls nX i=1 δy(i)1u(i)=l1x(i)∈Is where {Is}S s=1 is a...

  3. [3]

    We plan to investigate further its potential on large datasets in future work

    Assemble a product kernel with all these distributional encodings The advantage of this strategy is that interactions can be accounted for via step (b), but this requires a large training set and increases the problem dimension. We plan to investigate further its potential on large datasets in future work. 24 B Additional numerical experiments To compleme...

Pith tools

Reviewed August 7, 2026 · model on record in the stance chip above.