Sequential sparse Gaussian process quantile regression

# Sequential Sparse Gaussian Process Quantile Regression
Source: [https://arxiv.org/html/2606.31284](https://arxiv.org/html/2606.31284)
\\newsiamremark

remarkRemark\\newsiamremarkhypothesisHypothesis\\newsiamthmclaimClaim\\newsiamremarkfactFact\\headersSequential Sparse Gaussian Process Quantile RegressionH\. Nicolas and O\. Le Maître

Hugo NicolasInria, Center for Applied Mathematics, École polytechnique, Institut Polytechnique de Paris, Palaiseau, France \(\)\.Olivier Le MaîtreCNRS, Center for Applied Mathematics, École polytechnique, Institut Polytechnique de Paris, Palaiseau, France \(\)\.

###### Abstract

Quantile regression aims to estimate the conditional quantiles of a response variable from observed data\. In a Bayesian setting, Gaussian process quantile regression provides uncertainty quantification but faces significant computational challenges due to the nonconjugacy of the asymmetric Laplace likelihood and the cost of posterior inference\. We develop a sparse Gaussian process framework in which the quantile function is represented through a reduced set of inducing variables and posterior inference is performed using a Laplace approximation\. A decomposition of the predictive uncertainty into conditional\-prior and posterior\-induced variance components is then exploited to drive two complementary adaptive mechanisms: inducing\-input infilling and data acquisition\. These mechanisms are combined within a sequential algorithm that allocates computational effort toward the dominant source of predictive uncertainty and adaptively controls model complexity\. Numerical experiments on benchmark problems demonstrate the accuracy of the Laplace approximation, the benefits of variance\-based inducing\-input placement, and the effectiveness of the proposed sequential enrichment strategy compared with predefined data\-acquisition strategies\.

###### keywords:

Quantile Regression, Gaussian Processes, Laplace Approximation, Rejection Sampling, Uncertainty Quantification

\{MSCcodes\}

62G08, 62F15, 60G15, 68T05

## 1Introduction

Quantile regression aims to estimate the conditional quantiles of a response variable given a set of input variables\. It provides a more complete characterization of the conditional distribution than mean regression, and is particularly valuable in applications involving uncertainty quantification, risk assessment, reliability analysis, and decision\-making under uncertainty\.

Quantile estimation is classically based on the asymmetric loss function introduced by Koenker and Bassett\[[11](https://arxiv.org/html/2606.31284#bib.bib11)\], whose minimizer coincides with the desired conditional quantile\. Building on this formulation, a wide range of quantile regression functional forms have been proposed, including linear models\[[11](https://arxiv.org/html/2606.31284#bib.bib11),[33](https://arxiv.org/html/2606.31284#bib.bib33)\], spline\-based methods\[[12](https://arxiv.org/html/2606.31284#bib.bib12),[29](https://arxiv.org/html/2606.31284#bib.bib29)\], random forests\[[16](https://arxiv.org/html/2606.31284#bib.bib16)\], neural networks\[[5](https://arxiv.org/html/2606.31284#bib.bib5),[10](https://arxiv.org/html/2606.31284#bib.bib10)\], kernel methods\[[28](https://arxiv.org/html/2606.31284#bib.bib28)\], and Gaussian processes \(GPs\)\[[23](https://arxiv.org/html/2606.31284#bib.bib23),[22](https://arxiv.org/html/2606.31284#bib.bib22)\]\. In the present work, we focus on GP quantile regression\.

GPs provide a flexible nonparametric framework for modeling unknown functions while naturally quantifying predictive uncertainty\. Within a Bayesian formulation, uncertainty is represented through a posterior distribution over the quantile function and can be propagated to predictions\. A common choice for Bayesian quantile regression is to combine a GP prior with an asymmetric Laplace likelihood, for which maximum likelihood estimation recovers the frequentist quantile regression solution\. This combination offers a principled way to learn conditional quantiles together with their associated uncertainty\.

Two challenges limit the practical deployment of Bayesian quantile regression with GPs\. First, posterior inference is analytically intractable because the asymmetric Laplace likelihood is not conjugate to the GP prior\. Second, scalable sparse formulations require the construction of an efficient inducing representation, including both the number and the location of the inducing inputs\.

Several approximations were proposed to address this intractability\. Sampling\-based methods, such as Markov chain Monte Carlo \(MCMC\), are flexible for quantile regression tasks\[[33](https://arxiv.org/html/2606.31284#bib.bib33),[14](https://arxiv.org/html/2606.31284#bib.bib14),[13](https://arxiv.org/html/2606.31284#bib.bib13)\]but can become computationally demanding for large datasets and complex models\. Deterministic alternatives instead approximate the posterior directly, through expectation propagation \(EP\), variational inference, or the Laplace approximation\. In the context of GP quantile regression, EP was considered by Boukouvalas et al\.\[[3](https://arxiv.org/html/2606.31284#bib.bib3)\], while Abeywardana and Ramos\[[1](https://arxiv.org/html/2606.31284#bib.bib1)\]proposed a variational approximation based on the location\-scale mixture representation of the asymmetric Laplace distribution\.

To improve scalability, sparse GP models introduce a reduced set of latent inducing variables that act as a compact representation of the latent function\. Recent GP quantile regression frameworks, such as the one proposed by Picheny et al\.\[[21](https://arxiv.org/html/2606.31284#bib.bib21)\], have adopted sparse variational GP formulations\. Sparse variational GPs approximate the posterior distribution over the inducing variables by a Gaussian distribution whose mean and covariance are typically optimized jointly\[[27](https://arxiv.org/html/2606.31284#bib.bib27),[30](https://arxiv.org/html/2606.31284#bib.bib30),[8](https://arxiv.org/html/2606.31284#bib.bib8),[15](https://arxiv.org/html/2606.31284#bib.bib15)\]\. While effective, the numberMMof inducing points is typically specified*a priori*and the dimension of the optimization problem grows quadratically withMM\.

In this work, we develop a sequential sparse GP framework in which a decomposition of predictive uncertainty provides a unified basis for posterior inference, adaptive model enrichment, and data acquisition\. The proposed methodology addresses two key challenges in sparse GP quantile regression: the efficient approximation of the posterior distribution over the inducing variables and the adaptive control of model complexity\.

The first contribution is a sparse Bayesian quantile regression formulation based on a Laplace approximation of the posterior distribution over the inducing variables\. Rather than optimizing both the mean vector and covariance matrix of an approximate posterior distribution, as is done in sparse variational GPs, we recast inference as the determination of a maximum*a posteriori*estimate and the inversion of the associated Hessian matrix\. This yields a tractable optimization problem whose number of dominant optimization variables decreases from𝒪​\(M2\)\\mathcal\{O\}\(M^\{2\}\)for variational formulations to𝒪​\(M\)\\mathcal\{O\}\(M\)in the proposed approach\.

The second contribution is an adaptive inducing\-input infilling strategy, as the quality of the inducing\-point set is crucial for the performance of sparse GPs\[[4](https://arxiv.org/html/2606.31284#bib.bib4),[32](https://arxiv.org/html/2606.31284#bib.bib32),[17](https://arxiv.org/html/2606.31284#bib.bib17)\]\. The inducing inputs are treated as adaptive model components rather than fixed design variables\. Inspired by Burt et al\.\[[4](https://arxiv.org/html/2606.31284#bib.bib4)\]and Ober et al\.\[[20](https://arxiv.org/html/2606.31284#bib.bib20)\], new inducing inputs are introduced sequentially by maximizing the integrated reduction in conditional\-prior variance, thereby determining both their locations and their number\.

The third contribution is an adaptive data\-acquisition strategy based on the posterior uncertainty induced by the inducing variables\. Additional training data are acquired preferentially in regions where posterior uncertainty over the inducing variables contributes most strongly to predictive uncertainty\.

These two enrichment mechanisms act on distinct components of the predictive uncertainty: inducing\-point infilling reduces the conditional\-prior variance, whereas data acquisition reduces the posterior\-induced variance\. Finally, both enrichment mechanisms are combined into a unified sequential algorithm\. A variance\-based switching criterion determines whether computational effort should be allocated to improving the inducing representation or to acquiring additional training data\. The resulting procedure adaptively balances model complexity and data availability while targeting the dominant source of predictive uncertainty\.

The remainder of the paper is organized as follows\.[Section2](https://arxiv.org/html/2606.31284#S2)introduces the sparse GP quantile regression formulation and the Laplace approximation\.[Section3](https://arxiv.org/html/2606.31284#S3)presents the inducing\-input infilling strategy, the adaptive data\-acquisition procedure, and the sequential enrichment algorithm\. Numerical experiments are reported in[Section4](https://arxiv.org/html/2606.31284#S4), and technical derivations are collected in the appendices\.

## 2Bayesian quantile regression

### 2\.1Inference problem

Letf:𝒳×Ω→ℝf:\\mathcal\{X\}\\times\\varOmega\\to\\mathbb\{R\}be a measurable function, where the input space𝒳⊂ℝD\\mathcal\{X\}\\subset\\mathbb\{R\}^\{D\}is compact and\(Ω,ℱ,ℙ\)\(\\varOmega,\\mathcal\{F\},\\mathbb\{P\}\)denotes the underlying probability space\. Definey​\(𝐱\):=f​\(𝐱,⋅\)y\(\\mathbf\{x\}\):=f\(\\mathbf\{x\},\\cdot\)as the random output of the stochastic processffat a deterministic input𝐱∈𝒳\\mathbf\{x\}\\in\\mathcal\{X\}\. For a prescribed quantile levelτ∈\(0,1\)\\tau\\in\(0,1\), theτ\\tau\-quantile of the distribution ofy​\(𝐱\)y\(\\mathbf\{x\}\)is defined asqτ​\(𝐱\)=inf\{ν∈ℝ:ℙ​\(y​\(𝐱\)≤ν\)≥τ\}q\_\{\\tau\}\(\\mathbf\{x\}\)=\\inf\\\{\\nu\\in\\mathbb\{R\}:\\mathbb\{P\}\\,\(y\(\\mathbf\{x\}\)\\leq\\nu\)\\geq\\tau\\\}\.

Consider a training dataset𝒟=\{\(yn,𝐱n\)\}n=1N\\mathcal\{D\}=\\\{\(\\mathrm\{y\}\_\{n\},\\mathbf\{x\}\_\{n\}\)\\\}\_\{n=1\}^\{N\}ofNNindependent realizationsyn∈ℝ\\mathrm\{y\}\_\{n\}\\in\\mathbb\{R\}of random variablesy​\(𝐱n\)y\(\\mathbf\{x\}\_\{n\}\)and their associated input locations𝐱n∈𝒳\\mathbf\{x\}\_\{n\}\\in\\mathcal\{X\}\. Equivalently,yn=f​\(𝐱n,ωn\)\\mathrm\{y\}\_\{n\}=f\(\\mathbf\{x\}\_\{n\},\\omega\_\{n\}\)for independent drawsω1,…,ωN\\omega\_\{1\},\\ldots,\\omega\_\{N\}from\(Ω,ℱ,ℙ\)\(\\varOmega,\\mathcal\{F\},\\mathbb\{P\}\)\. Quantile regression aims to learn the unknownτ\\tau\-quantile functionqτ:𝒳→ℝq\_\{\\tau\}:\\mathcal\{X\}\\to\\mathbb\{R\}from the training data\.

In this work, we adopt a Bayesian perspective on the quantile regression problem\. Bayesian inference specifies a probabilistic model for the data\-generating process and offers a principled framework for updating the model with new observations, according to Bayes’ theorem\. The deviationε\\varepsilonfrom theτ\\tau\-quantile at an input𝐱∈𝒳\\mathbf\{x\}\\in\\mathcal\{X\}is a random variable modeled as

\(1\)ε​\(𝐱\)=y​\(𝐱\)−qτ​\(𝐱\)\.\\varepsilon\(\\mathbf\{x\}\)=y\(\\mathbf\{x\}\)\-q\_\{\\tau\}\(\\mathbf\{x\}\)\.
Following Yu and Moyeed\[[33](https://arxiv.org/html/2606.31284#bib.bib33)\], the deviation is assumed to follow an asymmetric Laplace distribution\. We consider the asymmetric Laplace density introduced by Yu and Zhang\[[34](https://arxiv.org/html/2606.31284#bib.bib34)\]and defined as

\(2\)fε​\(ν;α\)=τ​\(1−τ\)α​exp⁡\(−ρτ​\(ν\)α\),f\_\{\\varepsilon\}\(\\nu\\,;\\alpha\)=\\frac\{\\tau\(1\-\\tau\)\}\{\\alpha\}\\exp\\left\(\-\\frac\{\\rho\_\{\\tau\}\(\\nu\)\}\{\\alpha\}\\right\),whereα\>0\\alpha\>0is a scale parameter andρτ\\rho\_\{\\tau\}denotes the quantile check function\. It is given by

\(3\)ρτ​\(ν\)=\(τ−𝟙\{ν≤0\}\)​ν,\\rho\_\{\\tau\}\(\\nu\)=\\left\(\\tau\-\\mathds\{1\}\_\{\\\{\\nu\\leq 0\\\}\}\\right\)\\nu,where𝟙\{ν≤0\}\\mathds\{1\}\_\{\\\{\\nu\\leq 0\\\}\}denotes the indicator function that equals11ifν≤0\\nu\\leq 0and0otherwise\. The check function is nonnegative\. The scale parameterα\\alphaencodes the spread of the conditional distribution about itsτ\\tau\-quantile\. We assume that it is constant across the input space𝒳\\mathcal\{X\}\. This assumption can be relaxed by allowing the scale parameter to vary with the input, as done by Picheny et al\.\[[21](https://arxiv.org/html/2606.31284#bib.bib21)\], who placed a GP prior over it\. Under the conditional independence assumption, the observations𝐲=\[y1,…,yN\]⊤\\mathbf\{y\}=\[\\mathrm\{y\}\_\{1\},\\ldots,\\mathrm\{y\}\_\{N\}\]^\{\\top\}are independent given the input locations𝐗=\[𝐱1,…,𝐱N\]\\mathbf\{X\}=\[\\mathbf\{x\}\_\{1\},\\ldots,\\mathbf\{x\}\_\{N\}\]and theτ\\tau\-quantile functionqτq\_\{\\tau\}\. Consequently, the likelihoodp​\(𝐲∣𝐗,qτ\)p\(\\mathbf\{y\}\\mid\\mathbf\{X\},q\_\{\\tau\}\)of the observed data factorizes as

\(4\)p​\(𝐲∣𝐗,qτ\)=∏n=1Nfε​\(yn−qτ​\(𝐱n\);α\)=\(τ​\(1−τ\)α\)N​exp⁡\(−1α​∑n=1Nρτ​\(yn−qτ​\(𝐱n\)\)\)\.p\(\\mathbf\{y\}\\mid\\mathbf\{X\},q\_\{\\tau\}\)=\\prod\_\{n=1\}^\{N\}f\_\{\\varepsilon\}\(\\mathrm\{y\}\_\{n\}\-q\_\{\\tau\}\(\\mathbf\{x\}\_\{n\}\)\\,;\\alpha\)=\\left\(\\frac\{\\tau\(1\-\\tau\)\}\{\\alpha\}\\right\)^\{N\}\\exp\\left\(\-\\frac\{1\}\{\\alpha\}\\sum\_\{n=1\}^\{N\}\\rho\_\{\\tau\}\(\\mathrm\{y\}\_\{n\}\-q\_\{\\tau\}\(\\mathbf\{x\}\_\{n\}\)\)\\right\)\.
We model the latentτ\\tau\-quantile function as a realization of a GP\. A GP is a collection of random variables, any finite number of which have a joint normal distribution\[[24](https://arxiv.org/html/2606.31284#bib.bib24), Section 2\.2\]\. It defines a distribution over functions, fully specified by its mean and covariance functions\. The covariance function, also known as thekernel, encodes the assumptions about the smoothness and structure of the latent function\. Without loss of generality, we consider a zero\-mean GP\. The prior distribution over theτ\\tau\-quantile function is then

\(5\)π​\(qτ\)=𝒢​𝒫​\(0,κ​\(⋅,⋅;𝝍\)\),\\pi\(q\_\{\\tau\}\)=\\mathcal\{GP\}\\big\(0,\\kappa\(\\cdot,\\cdot\\,;\\boldsymbol\{\\psi\}\)\\big\),whereκ​\(⋅,⋅;𝝍\)\\kappa\(\\cdot,\\cdot\\,;\\boldsymbol\{\\psi\}\)denotes the kernel with hyperparameters𝝍∈Ψ\\boldsymbol\{\\psi\}\\in\\varPsi\. In the following, we adopt the anisotropic Matérn5/25/2kernel, defined as

\(6\)κ​\(𝐱,𝐱′;𝝍\)=σs2​\(1\+5​r​\(𝐱,𝐱′\)\+53​r2​\(𝐱,𝐱′\)\)​exp⁡\(−5​r​\(𝐱,𝐱′\)\),\\kappa\(\\mathbf\{x\},\\mathbf\{x\}^\{\\prime\}\\,;\\boldsymbol\{\\psi\}\)=\\sigma\_\{\\mathrm\{s\}\}^\{2\}\\left\(1\+\\sqrt\{5\}\\,r\(\\mathbf\{x\},\\mathbf\{x\}^\{\\prime\}\)\+\\frac\{5\}\{3\}r^\{2\}\(\\mathbf\{x\},\\mathbf\{x\}^\{\\prime\}\)\\right\)\\exp\\left\(\-\\sqrt\{5\}\\,r\(\\mathbf\{x\},\\mathbf\{x\}^\{\\prime\}\)\\right\),whereσs2∈ℝ\>0\\sigma\_\{\\text\{s\}\}^\{2\}\\in\\mathbb\{R\}\_\{\>0\}denotes the signal variance, and

\(7\)r​\(𝐱,𝐱′\)=\(𝐱−𝐱′\)⊤​𝚲−1​\(𝐱−𝐱′\)r\(\\mathbf\{x\},\\mathbf\{x\}^\{\\prime\}\)=\\sqrt\{\(\\mathbf\{x\}\-\\mathbf\{x\}^\{\\prime\}\)^\{\\top\}\\boldsymbol\{\\Lambda\}^\{\-1\}\(\\mathbf\{x\}\-\\mathbf\{x\}^\{\\prime\}\)\}is the anisotropic Euclidean distance, with𝚲:=diag⁡\(ℓ12,…,ℓD2\)\\boldsymbol\{\\Lambda\}:=\\operatorname\{diag\}\(\\ell\_\{1\}^\{2\},\\ldots,\\ell\_\{D\}^\{2\}\)a diagonal matrix of positive lengthscalesℓd∈ℝ\>0\\ell\_\{d\}\\in\\mathbb\{R\}\_\{\>0\}for each input dimensiond=1,…,Dd=1,\\dots,D\. The hyperparameter vector is thus given by𝝍=\[σs2,ℓ1,…,ℓD\]⊤\\boldsymbol\{\\psi\}=\[\\sigma\_\{\\text\{s\}\}^\{2\},\\ell\_\{1\},\\ldots,\\ell\_\{D\}\]^\{\\top\}\. In the remainder of this work, we simplify the notation by omitting hyperparameter dependencies that are not essential to the discussion\.

Bayes’ theorem defines the posterior distribution over theτ\\tau\-quantile function conditional on the training data𝒟=\(𝐲,𝐗\)\\mathcal\{D\}=\(\\mathbf\{y\},\\mathbf\{X\}\)as

\(8\)p​\(qτ∣𝐲,𝐗\)=p​\(𝐲∣𝐗,qτ\)​π​\(qτ\)p​\(𝐲∣𝐗\)\.p\(q\_\{\\tau\}\\mid\\mathbf\{y\},\\mathbf\{X\}\)=\\frac\{p\(\\mathbf\{y\}\\mid\\mathbf\{X\},q\_\{\\tau\}\)\\,\\pi\(q\_\{\\tau\}\)\}\{p\(\\mathbf\{y\}\\mid\\mathbf\{X\}\)\}\.The marginal likelihoodp​\(𝐲∣𝐗\)p\(\\mathbf\{y\}\\mid\\mathbf\{X\}\)ensures the normalization of the posterior\. It is defined as the expectation of the likelihoodp​\(𝐲∣𝐗,qτ\)p\(\\mathbf\{y\}\\mid\\mathbf\{X\},q\_\{\\tau\}\)under the GP priorπ​\(qτ\)\\pi\(q\_\{\\tau\}\)\.

### 2\.2Sparse formulation

We introduce latent auxiliary variables that act as anchors for the regression\. Each auxiliary variable represents an evaluation of the latentτ\\tau\-quantile functionqτq\_\{\\tau\}at a location𝐳m∈𝒳\\mathbf\{z\}\_\{m\}\\in\\mathcal\{X\}\. Denoting the set of these locations by𝒵=\{𝐳m\}m=1M\\mathcal\{Z\}=\\\{\\mathbf\{z\}\_\{m\}\\\}\_\{m=1\}^\{M\}, we define these auxiliary variables as𝒖=\[qτ​\(𝐳1\),…,qτ​\(𝐳M\)\]⊤\\boldsymbol\{u\}=\[q\_\{\\tau\}\(\\mathbf\{z\}\_\{1\}\),\\ldots,q\_\{\\tau\}\(\\mathbf\{z\}\_\{M\}\)\]^\{\\top\}\. In the remainder, we refer to the elements of𝒖\\boldsymbol\{u\}asinducing variables, and the locations in𝒵\\mathcal\{Z\}asinducing inputs\. Crucially, their number satisfiesM≪NM\\ll N, which ensures scalable inference for large training datasets\. They may be taken either as a subset of the training inputs𝐗\\mathbf\{X\}or as an entirely distinct set of locations\. Their selection is a key aspect of the proposed method, as we shall see in[Section3\.1](https://arxiv.org/html/2606.31284#S3.SS1)\. Until then, we assume that the inducing input locations in𝒵\\mathcal\{Z\}are fixed and known\.

#### 2\.2\.1Prior distribution

Denote by𝐮\\mathbf\{u\}a realization of the random vector𝒖\\boldsymbol\{u\}\. Since theτ\\tau\-quantile function carries a GP prior, the inducing variables, conditional on the inducing inputs𝐙=\[𝐳1,…,𝐳M\]\\mathbf\{Z\}=\[\\mathbf\{z\}\_\{1\},\\ldots,\\mathbf\{z\}\_\{M\}\], follow a multivariate normal distribution\. Its probability density function \(PDF\) is

\(9\)π​\(𝐮∣𝐙\)=𝒩​\(𝐮∣𝟎,K​\(𝐙,𝐙\)\)=1\(2​π\)M2​\|K​\(𝐙,𝐙\)\|12​exp⁡\(−12​𝐮⊤​K​\(𝐙,𝐙\)−1​𝐮\),\\pi\(\\mathbf\{u\}\\mid\\mathbf\{Z\}\)=\\mathcal\{N\}\\big\(\\mathbf\{u\}\\mid\\boldsymbol\{0\},K\(\\mathbf\{Z\},\\mathbf\{Z\}\)\\big\)=\\frac\{1\}\{\(2\\pi\)^\{\\frac\{M\}\{2\}\}\\left\\lvert K\(\\mathbf\{Z\},\\mathbf\{Z\}\)\\right\\rvert^\{\\frac\{1\}\{2\}\}\}\\exp\\left\(\-\\frac\{1\}\{2\}\\mathbf\{u\}^\{\\top\}K\(\\mathbf\{Z\},\\mathbf\{Z\}\)^\{\-1\}\\mathbf\{u\}\\right\),where\|⋅\|\\left\\lvert\\cdot\\right\\rvertdenotes the determinant operator, andK​\(𝐙,𝐙\)K\(\\mathbf\{Z\},\\mathbf\{Z\}\)is the covariance matrix evaluated at the inducing inputs, with entries\[K​\(𝐙,𝐙\)\]i​j=κ​\(𝐳i,𝐳j;𝝍\)\[K\(\\mathbf\{Z\},\\mathbf\{Z\}\)\]\_\{ij\}=\\kappa\(\\mathbf\{z\}\_\{i\},\\mathbf\{z\}\_\{j\}\\,;\\boldsymbol\{\\psi\}\)fori,j=1,…,Mi,j=1,\\ldots,M\. The prior over theτ\\tau\-quantile function, conditional on the inducing variables, is

\(10\)π​\(qτ∣𝐮,𝐙\)=𝒢​𝒫​\(μ​\(⋅;𝐮,𝐙\),Σz2​\(⋅,⋅;𝐙\)\)\.\\pi\(q\_\{\\tau\}\\mid\\mathbf\{u\},\\mathbf\{Z\}\)=\\mathcal\{GP\}\\big\(\\mu\(\\cdot\\,;\\mathbf\{u\},\\mathbf\{Z\}\),\\Sigma\_\{z\}^\{2\}\(\\cdot,\\cdot\\,;\\mathbf\{Z\}\)\\big\)\.The mean and covariance functions are, respectively, given by

\(11\)μ​\(𝐱;𝐮,𝐙\)\\displaystyle\\mu\(\\mathbf\{x\}\\,;\\mathbf\{u\},\\mathbf\{Z\}\)=K​\(𝐱,𝐙\)​K​\(𝐙,𝐙\)−1​𝐮,\\displaystyle=K\(\\mathbf\{x\},\\mathbf\{Z\}\)K\(\\mathbf\{Z\},\\mathbf\{Z\}\)^\{\-1\}\\mathbf\{u\},\(12\)Σz2​\(𝐱,𝐱′;𝐙\)\\displaystyle\\Sigma\_\{z\}^\{2\}\(\\mathbf\{x\},\\mathbf\{x\}^\{\\prime\}\\,;\\mathbf\{Z\}\)=K​\(𝐱,𝐱′\)−K​\(𝐱,𝐙\)​K​\(𝐙,𝐙\)−1​K​\(𝐙,𝐱′\)\.\\displaystyle=K\(\\mathbf\{x\},\\mathbf\{x\}^\{\\prime\}\)\-K\(\\mathbf\{x\},\\mathbf\{Z\}\)K\(\\mathbf\{Z\},\\mathbf\{Z\}\)^\{\-1\}K\(\\mathbf\{Z\},\\mathbf\{x\}^\{\\prime\}\)\.

#### 2\.2\.2Posterior distribution

We specify the joint posterior over theτ\\tau\-quantile function and the inducing variables as

\(13\)p​\(qτ,𝐮∣𝐙,𝐲,𝐗\)=π​\(qτ∣𝐮,𝐙\)​p​\(𝐮∣𝐙,𝐲,𝐗\)\.p\(q\_\{\\tau\},\\mathbf\{u\}\\mid\\mathbf\{Z\},\\mathbf\{y\},\\mathbf\{X\}\)=\\pi\(q\_\{\\tau\}\\mid\\mathbf\{u\},\\mathbf\{Z\}\)\\,p\(\\mathbf\{u\}\\mid\\mathbf\{Z\},\\mathbf\{y\},\\mathbf\{X\}\)\.The first term on the right\-hand side is the conditional prior over theτ\\tau\-quantile function, given in[Eq\.10](https://arxiv.org/html/2606.31284#S2.E10)\. The second term is the posterior over the inducing variables, conditional on the training data\. By Bayes’ theorem, it decomposes as

\(14\)p​\(𝐮∣𝐙,𝐲,𝐗\)=p​\(𝐲∣𝐗,𝐮,𝐙\)​π​\(𝐮∣𝐙\)p​\(𝐲∣𝐗\)\.p\(\\mathbf\{u\}\\mid\\mathbf\{Z\},\\mathbf\{y\},\\mathbf\{X\}\)=\\frac\{p\(\\mathbf\{y\}\\mid\\mathbf\{X\},\\mathbf\{u\},\\mathbf\{Z\}\)\\,\\pi\(\\mathbf\{u\}\\mid\\mathbf\{Z\}\)\}\{p\(\\mathbf\{y\}\\mid\\mathbf\{X\}\)\}\.The priorπ​\(𝐮∣𝐙\)\\pi\(\\mathbf\{u\}\\mid\\mathbf\{Z\}\)over the inducing variables is given in[Eq\.9](https://arxiv.org/html/2606.31284#S2.E9)\. The conditional likelihood

\(15\)p​\(𝐲∣𝐗,𝐮,𝐙\)=𝔼π​\(qτ∣𝐮,𝐙\)​\[p​\(𝐲∣𝐗,qτ\)\],p\(\\mathbf\{y\}\\mid\\mathbf\{X\},\\mathbf\{u\},\\mathbf\{Z\}\)=\\mathbb\{E\}\_\{\\pi\(q\_\{\\tau\}\\mid\\mathbf\{u\},\\mathbf\{Z\}\)\}\\big\[p\(\\mathbf\{y\}\\mid\\mathbf\{X\},q\_\{\\tau\}\)\\big\],marginalizes the data likelihoodp​\(𝐲∣𝐗,qτ\)p\(\\mathbf\{y\}\\mid\\mathbf\{X\},q\_\{\\tau\}\), defined in[Eq\.4](https://arxiv.org/html/2606.31284#S2.E4), over the conditional prior \([10](https://arxiv.org/html/2606.31284#S2.E10)\)\.

The posterior quantile modelp​\(qτ∣𝐲,𝐗,𝐙\)p\(q\_\{\\tau\}\\mid\\mathbf\{y\},\\mathbf\{X\},\\mathbf\{Z\}\)can be obtained by marginalizing the joint posteriorp​\(qτ,𝐮∣𝐙,𝐲,𝐗\)p\(q\_\{\\tau\},\\mathbf\{u\}\\mid\\mathbf\{Z\},\\mathbf\{y\},\\mathbf\{X\}\)over the posteriorp​\(𝐮∣𝐙,𝐲,𝐗\)p\(\\mathbf\{u\}\\mid\\mathbf\{Z\},\\mathbf\{y\},\\mathbf\{X\}\)\. However, in a significant departure from classical GP regression, the training values consist of realizations of the response variabley​\(𝐱\)y\(\\mathbf\{x\}\)rather than direct observations of the conditionalτ\\tau\-quantile\. Formally, the asymmetric Laplace likelihood is not conjugate to the GP prior, rendering the posteriorp​\(𝐮∣𝐙,𝐲,𝐗\)p\(\\mathbf\{u\}\\mid\\mathbf\{Z\},\\mathbf\{y\},\\mathbf\{X\}\)over the inducing variables analytically intractable\. Approximating this posterior restores analytical tractability and is the subject of[Section2\.3](https://arxiv.org/html/2606.31284#S2.SS3)\.

### 2\.3Approximate inference

#### 2\.3\.1Laplace approximation

We approximate the true posteriorp​\(𝐮∣𝐙,𝐲,𝐗\)p\(\\mathbf\{u\}\\mid\\mathbf\{Z\},\\mathbf\{y\},\\mathbf\{X\}\)over the inducing variables using the Laplace approximation, a technique that replaces an intractable posterior with a Gaussian centered at its MAP estimate\[[18](https://arxiv.org/html/2606.31284#bib.bib18), Section 4\.6\.8\]\.

Given the training dataset𝒟=\(𝐲,𝐗\)\\mathcal\{D\}=\(\\mathbf\{y\},\\mathbf\{X\}\), the MAP estimate𝐮^\\hat\{\\mathbf\{u\}\}of the inducing variables can be found by maximizing the logarithm of the posteriorp​\(𝐮∣𝐙,𝐲,𝐗\)p\(\\mathbf\{u\}\\mid\\mathbf\{Z\},\\mathbf\{y\},\\mathbf\{X\}\)\. This is equivalent to maximizing the logarithm of the joint densityp​\(𝐮,𝐲∣𝐗,𝐙\)p\(\\mathbf\{u\},\\mathbf\{y\}\\mid\\mathbf\{X\},\\mathbf\{Z\}\), such that

\(16\)𝐮^=arg⁡max𝐮∈ℝM⁡𝒥N​\(𝐮\),with𝒥N​\(𝐮\):=log⁡p​\(𝐮,𝐲∣𝐗,𝐙\)=log⁡p​\(𝐲∣𝐗,𝐮,𝐙\)\+log⁡π​\(𝐮∣𝐙\)\.\\hat\{\\mathbf\{u\}\}=\\arg\\max\_\{\\mathbf\{u\}\\in\\mathbb\{R\}^\{M\}\}\\mathcal\{J\}\_\{N\}\(\\mathbf\{u\}\),\\quad\\text\{with\}\\quad\\mathcal\{J\}\_\{N\}\(\\mathbf\{u\}\):=\\log p\(\\mathbf\{u\},\\mathbf\{y\}\\mid\\mathbf\{X\},\\mathbf\{Z\}\)=\\log p\(\\mathbf\{y\}\\mid\\mathbf\{X\},\\mathbf\{u\},\\mathbf\{Z\}\)\+\\log\\pi\(\\mathbf\{u\}\\mid\\mathbf\{Z\}\)\.
The conditional likelihoodp​\(𝐲∣𝐗,𝐮,𝐙\)p\(\\mathbf\{y\}\\mid\\mathbf\{X\},\\mathbf\{u\},\\mathbf\{Z\}\), given in[Eq\.15](https://arxiv.org/html/2606.31284#S2.E15), does not admit a closed\-form expression and would therefore require approximation by sampling\. Repeated sampling would be prohibitively expensive given the sequential refinement strategy proposed in the present work\. To overcome this limitation, we propose to use a surrogate𝒥~N\\widetilde\{\\mathcal\{J\}\}\_\{N\}of the log\-joint density𝒥N​\(𝐮\)=log⁡\(𝔼π​\(qτ∣𝐮,𝐙\)​\[p​\(𝐲∣𝐗,qτ\)\]\)\+log⁡π​\(𝐮∣𝐙\)\\mathcal\{J\}\_\{N\}\(\\mathbf\{u\}\)=\\log\\left\(\\mathbb\{E\}\_\{\\pi\(q\_\{\\tau\}\\mid\\mathbf\{u\},\\mathbf\{Z\}\)\}\\big\[p\(\\mathbf\{y\}\\mid\\mathbf\{X\},q\_\{\\tau\}\)\\big\]\\right\)\+\\log\\pi\(\\mathbf\{u\}\\mid\\mathbf\{Z\}\)\. The surrogate is obtained by interchanging the expectation and the logarithm, yielding

\(17\)𝒥~N​\(𝐮\):=𝔼π​\(qτ∣𝐮,𝐙\)​\[log⁡p​\(𝐲∣𝐗,qτ\)\]\+log⁡π​\(𝐮∣𝐙\)\.\\widetilde\{\\mathcal\{J\}\}\_\{N\}\(\\mathbf\{u\}\):=\\mathbb\{E\}\_\{\\pi\(q\_\{\\tau\}\\mid\\mathbf\{u\},\\mathbf\{Z\}\)\}\\big\[\\log p\(\\mathbf\{y\}\\mid\\mathbf\{X\},q\_\{\\tau\}\)\\big\]\+\\log\\pi\(\\mathbf\{u\}\\mid\\mathbf\{Z\}\)\.By Jensen’s inequality,

\(18\)𝔼π​\(qτ∣𝐮,𝐙\)​\[log⁡p​\(𝐲∣𝐗,qτ\)\]≤log⁡\(𝔼π​\(qτ∣𝐮,𝐙\)​\[p​\(𝐲∣𝐗,qτ\)\]\),\\mathbb\{E\}\_\{\\pi\(q\_\{\\tau\}\\mid\\mathbf\{u\},\\mathbf\{Z\}\)\}\\big\[\\log p\(\\mathbf\{y\}\\mid\\mathbf\{X\},q\_\{\\tau\}\)\\big\]\\leq\\log\\left\(\\mathbb\{E\}\_\{\\pi\(q\_\{\\tau\}\\mid\\mathbf\{u\},\\mathbf\{Z\}\)\}\\big\[p\(\\mathbf\{y\}\\mid\\mathbf\{X\},q\_\{\\tau\}\)\\big\]\\right\),so this surrogate is not arbitrary: it is a lower bound on the exact log\-joint objective and corresponds to maximizing a tractable, conservative approximation of it\.[AppendixA](https://arxiv.org/html/2606.31284#A1)shows that the normalized gap between the two objectives is controlled by the average conditional\-prior variance at the training inputs, thereby connecting the accuracy of the surrogate to the inducing\-input infilling strategy\. Although the effect of the surrogate on the location of the MAP estimate𝐮^\\hat\{\\mathbf\{u\}\}is not analyzed theoretically in this work, the numerical results support the practical effectiveness of the surrogate\. In addition, the surrogate has a closed\-form expression that enables the derivation of its gradient and Hessian, provided in[AppendixB](https://arxiv.org/html/2606.31284#A2)\. In fact, the Gaussian convolution of the check functionρτ\\rho\_\{\\tau\}smooths out its nondifferentiability, which, combined with theC∞C^\{\\infty\}smoothness of the Gaussian prior in𝐮\\mathbf\{u\}and the linearity in𝐮\\mathbf\{u\}of the conditional\-prior mean functionμ\\mu, renders the approximate log\-joint density𝒥~N​\(𝐮\)\\widetilde\{\\mathcal\{J\}\}\_\{N\}\(\\mathbf\{u\}\)of classC∞C^\{\\infty\}\. These characteristics enable solving[Eq\.16](https://arxiv.org/html/2606.31284#S2.E16), with𝒥~N\\widetilde\{\\mathcal\{J\}\}\_\{N\}in place of𝒥N\\mathcal\{J\}\_\{N\}, using efficient deterministic gradient\-based optimization algorithms\. The availability of the Hessian can be further exploited to accelerate the convergence of the optimization\. For these reasons, the approximate log\-density𝒥~N\\widetilde\{\\mathcal\{J\}\}\_\{N\}is used in the rest of the derivation\.

A second\-order Taylor expansion of the log\-posterior about𝐮^\\hat\{\\mathbf\{u\}\}yields the Laplace approximation

\(19\)p​\(𝐮∣𝐙,𝐲,𝐗\)≃p^​\(𝐮∣𝐙,𝐲,𝐗\):=𝒩​\(𝐮∣𝐮^,𝐂^\),p\(\\mathbf\{u\}\\mid\\mathbf\{Z\},\\mathbf\{y\},\\mathbf\{X\}\)\\simeq\\hat\{p\}\(\\mathbf\{u\}\\mid\\mathbf\{Z\},\\mathbf\{y\},\\mathbf\{X\}\):=\\mathcal\{N\}\\big\(\\mathbf\{u\}\\mid\\hat\{\\mathbf\{u\}\},\\hat\{\\mathbf\{C\}\}\\big\),where𝐂^\\hat\{\\mathbf\{C\}\}denotes the inverse of the negative log\-posterior Hessian evaluated at the MAP estimate:

\(20\)𝐂^=\[−𝐇N​\(𝐮^\)\]−1,with𝐇N​\(𝐮^\):=∇𝐮2𝒥~N​\(𝐮\)\|𝐮=𝐮^\.\\hat\{\\mathbf\{C\}\}=\[\-\\mathbf\{H\}\_\{N\}\(\\hat\{\\mathbf\{u\}\}\)\]^\{\-1\},\\quad\\text\{with\}\\quad\\mathbf\{H\}\_\{N\}\(\\hat\{\\mathbf\{u\}\}\):=\\nabla\_\{\\mathbf\{u\}\}^\{2\}\\,\\widetilde\{\\mathcal\{J\}\}\_\{N\}\(\\mathbf\{u\}\)\|\_\{\\mathbf\{u\}=\\hat\{\\mathbf\{u\}\}\}\.
In practice, the GP hyperparameters𝝍\\boldsymbol\{\\psi\}and the asymmetric Laplace scale parameterα\\alphaenter the MAP objective\. We therefore reintroduce them explicitly into the notation and optimize them alongside𝐮^\\hat\{\\mathbf\{u\}\}, leading to the following bilevel optimization problem:

\(21a\)\(𝝍^,α^\)\\displaystyle\\big\(\\hat\{\\boldsymbol\{\\psi\}\},\\hat\{\\alpha\}\\big\)=arg⁡max\(𝝍,α\)∈Ψ×ℝ\>0ℱN​\(𝝍,α\),\\displaystyle=\\mathop\{\\arg\\max\}\_\{\(\\boldsymbol\{\\psi\},\\alpha\)\\in\\varPsi\\times\\mathbb\{R\}\_\{\>0\}\}\\mathcal\{F\}\_\{N\}\(\\boldsymbol\{\\psi\},\\alpha\),\(21b\)subject to𝐮^​\(𝝍,α\)\\displaystyle\\text\{subject to\}\\quad\\hat\{\\mathbf\{u\}\}\(\\boldsymbol\{\\psi\},\\alpha\)=arg⁡max𝐮∈ℝM⁡𝒥~N​\(𝐮;𝝍,α\),\\displaystyle=\\arg\\max\_\{\\mathbf\{u\}\\in\\mathbb\{R\}^\{M\}\}\\widetilde\{\\mathcal\{J\}\}\_\{N\}\(\\mathbf\{u\}\\,;\\boldsymbol\{\\psi\},\\alpha\),whereℱN\\mathcal\{F\}\_\{N\}denotes the log\-marginal likelihood, expressed as

\(22a\)ℱN​\(𝝍,α\)\\displaystyle\\mathcal\{F\}\_\{N\}\(\\boldsymbol\{\\psi\},\\alpha\):=log⁡p​\(𝐲∣𝐗,𝝍,α\)\\displaystyle:=\\log p\(\\mathbf\{y\}\\mid\\mathbf\{X\},\\boldsymbol\{\\psi\},\\alpha\)\(22b\)=log​∫ℝMexp⁡\(𝒥N​\(𝐮;𝝍,α\)\)​𝑑𝐮\\displaystyle=\\log\\int\_\{\\mathbb\{R\}^\{M\}\}\\exp\\left\(\\mathcal\{J\}\_\{N\}\(\\mathbf\{u\}\\,;\\boldsymbol\{\\psi\},\\alpha\)\\right\)\\,d\\mathbf\{u\}\(22c\)≃𝒥~N​\(𝐮^​\(𝝍,α\);𝝍,α\)−12​log⁡\|−𝐇N​\(𝐮^​\(𝝍,α\);𝝍,α\)\|\+Cst,\\displaystyle\\simeq\\widetilde\{\\mathcal\{J\}\}\_\{N\}\(\\hat\{\\mathbf\{u\}\}\(\\boldsymbol\{\\psi\},\\alpha\)\\,;\\boldsymbol\{\\psi\},\\alpha\)\-\\frac\{1\}\{2\}\\log\\left\\lvert\-\\mathbf\{H\}\_\{N\}\(\\hat\{\\mathbf\{u\}\}\(\\boldsymbol\{\\psi\},\\alpha\)\\,;\\boldsymbol\{\\psi\},\\alpha\)\\right\\rvert\+\\mathrm\{C^\{st\}\},with the constant termCst\\mathrm\{C^\{st\}\}independent of𝝍\\boldsymbol\{\\psi\}andα\\alpha\. To improve numerical conditioning, the inner\-problem optimization[Eq\.21b](https://arxiv.org/html/2606.31284#S2.E21.2)is carried out in the whitened parameterization𝐮=𝐋𝐯\\mathbf\{u\}=\\mathbf\{L\}\\mathbf\{v\}\[[19](https://arxiv.org/html/2606.31284#bib.bib19),[9](https://arxiv.org/html/2606.31284#bib.bib9)\], where𝐋\\mathbf\{L\}is the lower Cholesky factor ofK​\(𝐙,𝐙\)K\(\\mathbf\{Z\},\\mathbf\{Z\}\), so that𝐯∼𝒩​\(𝟎,𝐈M\)\\mathbf\{v\}\\sim\\mathcal\{N\}\(\\mathbf\{0\},\\mathbf\{I\}\_\{M\}\)under the prior\.

#### 2\.3\.2Predictive quantile distribution

The MAP estimate𝐮^\\hat\{\\mathbf\{u\}\}and the covariance matrix𝐂^\\hat\{\\mathbf\{C\}\}in the Laplace approximation are obtained for a specific training dataset𝒟=\(𝐲,𝐗\)\\mathcal\{D\}=\(\\mathbf\{y\},\\mathbf\{X\}\)\. We make this dependence explicit by writing these quantities as𝐮^​\(𝐲,𝐗\)\\hat\{\\mathbf\{u\}\}\(\\mathbf\{y\},\\mathbf\{X\}\)and𝐂^​\(𝐲,𝐗\)\\hat\{\\mathbf\{C\}\}\(\\mathbf\{y\},\\mathbf\{X\}\), respectively\. Substituting the Laplace approximationp^​\(𝐮∣𝐙,𝐲,𝐗\)\\hat\{p\}\(\\mathbf\{u\}\\mid\\mathbf\{Z\},\\mathbf\{y\},\\mathbf\{X\}\)for the true posteriorp​\(𝐮∣𝐙,𝐲,𝐗\)p\(\\mathbf\{u\}\\mid\\mathbf\{Z\},\\mathbf\{y\},\\mathbf\{X\}\)yields the following approximate posterior predictive distribution over theτ\\tau\-quantile function:

\(23\)p^​\(qτ∣𝐲,𝐗,𝐙\)=𝔼p^​\(𝐮∣𝐙,𝐲,𝐗\)​\[π​\(qτ∣𝐮,𝐙\)\]=𝒢​𝒫​\(q^τ​\(⋅;𝐲,𝐗,𝐙\),Σz2​\(⋅,⋅;𝐙\)\+Σu2​\(⋅,⋅;𝐲,𝐗,𝐙\)\),\\hat\{p\}\(q\_\{\\tau\}\\mid\\mathbf\{y\},\\mathbf\{X\},\\mathbf\{Z\}\)=\\mathbb\{E\}\_\{\\hat\{p\}\(\\mathbf\{u\}\\mid\\mathbf\{Z\},\\mathbf\{y\},\\mathbf\{X\}\)\}\\big\[\\pi\(q\_\{\\tau\}\\mid\\mathbf\{u\},\\mathbf\{Z\}\)\\big\]=\\mathcal\{GP\}\\big\(\\hat\{q\}\_\{\\tau\}\(\\cdot\\,;\\mathbf\{y\},\\mathbf\{X\},\\mathbf\{Z\}\),\\Sigma\_\{z\}^\{2\}\(\\cdot,\\cdot\\,;\\mathbf\{Z\}\)\+\\Sigma\_\{u\}^\{2\}\(\\cdot,\\cdot\\,;\\mathbf\{y\},\\mathbf\{X\},\\mathbf\{Z\}\)\\big\),whereq^τ\\hat\{q\}\_\{\\tau\}denotes the posterior predictive mean function, and the posterior predictive covariance function decomposes into the conditional\-prior covarianceΣz2\\Sigma\_\{z\}^\{2\}and the posterior\-induced covarianceΣu2\\Sigma\_\{u\}^\{2\}\. Their explicit expressions are given by

\(24\)q^τ​\(𝐱;𝐲,𝐗,𝐙\)\\displaystyle\\hat\{q\}\_\{\\tau\}\(\\mathbf\{x\}\\,;\\mathbf\{y\},\\mathbf\{X\},\\mathbf\{Z\}\)=K​\(𝐱,𝐙\)​K​\(𝐙,𝐙\)−1​𝐮^​\(𝐲,𝐗\),\\displaystyle=K\(\\mathbf\{x\},\\mathbf\{Z\}\)K\(\\mathbf\{Z\},\\mathbf\{Z\}\)^\{\-1\}\\hat\{\\mathbf\{u\}\}\(\\mathbf\{y\},\\mathbf\{X\}\),\(25\)Σz2​\(𝐱,𝐱′;𝐙\)\\displaystyle\\Sigma\_\{z\}^\{2\}\(\\mathbf\{x\},\\mathbf\{x\}^\{\\prime\}\\,;\\mathbf\{Z\}\)=K​\(𝐱,𝐱′\)−K​\(𝐱,𝐙\)​K​\(𝐙,𝐙\)−1​K​\(𝐙,𝐱′\),\\displaystyle=K\(\\mathbf\{x\},\\mathbf\{x\}^\{\\prime\}\)\-K\(\\mathbf\{x\},\\mathbf\{Z\}\)K\(\\mathbf\{Z\},\\mathbf\{Z\}\)^\{\-1\}K\(\\mathbf\{Z\},\\mathbf\{x\}^\{\\prime\}\),\(26\)Σu2​\(𝐱,𝐱′;𝐲,𝐗,𝐙\)\\displaystyle\\Sigma^\{2\}\_\{u\}\(\\mathbf\{x\},\\mathbf\{x\}^\{\\prime\}\\,;\\mathbf\{y\},\\mathbf\{X\},\\mathbf\{Z\}\)=K​\(𝐱,𝐙\)​K​\(𝐙,𝐙\)−1​𝐂^​\(𝐲,𝐗\)​K​\(𝐙,𝐙\)−1​K​\(𝐙,𝐱′\)\.\\displaystyle=K\(\\mathbf\{x\},\\mathbf\{Z\}\)K\(\\mathbf\{Z\},\\mathbf\{Z\}\)^\{\-1\}\\hat\{\\mathbf\{C\}\}\(\\mathbf\{y\},\\mathbf\{X\}\)K\(\\mathbf\{Z\},\\mathbf\{Z\}\)^\{\-1\}K\(\\mathbf\{Z\},\\mathbf\{x\}^\{\\prime\}\)\.

## 3Sequential enrichment

We adopt the integrated mean squared error \(IMSE\)\[[26](https://arxiv.org/html/2606.31284#bib.bib26)\]as a performance metric to quantify predictive accuracy\. For the present quantile model, it is expressed as

\(27a\)IMSE=∫𝒳𝔼p^​\(qτ∣𝐲,𝐗,𝐙\)​\[\(qτ​\(𝐱\)−qτ,true​\(𝐱\)\)2\]​𝑑𝐱\\displaystyle=\\int\_\{\\mathcal\{X\}\}\\mathbb\{E\}\_\{\\hat\{p\}\(q\_\{\\tau\}\\mid\\mathbf\{y\},\\mathbf\{X\},\\mathbf\{Z\}\)\}\\big\[\\left\(q\_\{\\tau\}\(\\mathbf\{x\}\)\-q\_\{\\tau,\\text\{true\}\}\(\\mathbf\{x\}\)\\right\)^\{2\}\\big\]\\,d\\mathbf\{x\}\(27b\)=∫𝒳\{Σz2​\(𝐱,𝐱;𝐙\)\+Σu2​\(𝐱,𝐱;𝐲,𝐗,𝐙\)⏟variance\+\(q^τ​\(𝐱;𝐲,𝐗,𝐙\)−qτ,true​\(𝐱\)\)2⏟best\-prediction squared error\}​𝑑𝐱,\\displaystyle=\\int\_\{\\mathcal\{X\}\}\\Big\\\{\\underbrace\{\\Sigma\_\{z\}^\{2\}\(\\mathbf\{x\},\\mathbf\{x\}\\,;\\mathbf\{Z\}\)\+\\Sigma\_\{u\}^\{2\}\(\\mathbf\{x\},\\mathbf\{x\}\\,;\\mathbf\{y\},\\mathbf\{X\},\\mathbf\{Z\}\)\}\_\{\\text\{variance\}\}\+\\underbrace\{\\left\(\\hat\{q\}\_\{\\tau\}\(\\mathbf\{x\}\\,;\\mathbf\{y\},\\mathbf\{X\},\\mathbf\{Z\}\)\-q\_\{\\tau,\\text\{true\}\}\(\\mathbf\{x\}\)\\right\)^\{2\}\}\_\{\\text\{best\-prediction squared error\}\}\\Big\\\}\\,d\\mathbf\{x\},whereqτ,true​\(𝐱\)q\_\{\\tau,\\text\{true\}\}\(\\mathbf\{x\}\)denotes the trueτ\\tau\-quantile of the distribution ofy​\(𝐱\)y\(\\mathbf\{x\}\)\.

The first variance term in[Eq\.27b](https://arxiv.org/html/2606.31284#S3.E27.2)is the conditional\-prior variance evaluated at input𝐱∈𝒳\\mathbf\{x\}\\in\\mathcal\{X\}\. With fixed GP hyperparameters, it is reduced by adding distinct inducing inputs\. Judicious inducing\-input placement accelerates this reduction\. The second variance term in[Eq\.27b](https://arxiv.org/html/2606.31284#S3.E27.2)reflects the variability in the quantile predictions induced by the posterior uncertainty over the inducing variables𝒖\\boldsymbol\{u\}\. This component is reduced by improving the estimation of the inducing variables\. With fixed hyperparameters, this can be achieved by augmenting the training dataset𝒟=\(𝐲,𝐗\)\\mathcal\{D\}=\(\\mathbf\{y\},\\mathbf\{X\}\)with additional observations\. Both variance terms can be computed in closed form, respectively through the covariance functions \([25](https://arxiv.org/html/2606.31284#S2.E25)\) and \([26](https://arxiv.org/html/2606.31284#S2.E26)\)\. The final term in[Eq\.27b](https://arxiv.org/html/2606.31284#S3.E27.2), referred to as the best\-prediction squared error, measures the squared deviation between the posterior predictive mean and the trueτ\\tau\-quantile\. It is typically not available since the trueτ\\tau\-quantile is unknown\.

We present a sequential enrichment strategy to improve the global accuracy of the posterior quantile model\. The strategy comprises two adaptive components, each targeting a distinct variance contribution in the IMSE \([27b](https://arxiv.org/html/2606.31284#S3.E27.2)\)\. The first component is an inducing\-input infilling procedure that maximizes the integrated reduction in conditional\-prior variance\. This yields a systematic and optimal placement of inducing inputs𝐙\\mathbf\{Z\}for a given training dataset𝒟\\mathcal\{D\}\. The second component concerns the acquisition of new observations of the stochastic processy​\(⋅\)y\(\\cdot\)\. It determines where training data should be acquired, based on the posterior\-induced variance of the posterior quantile model, whose complexity is governed by the number and placement of the inducing inputs𝐙\\mathbf\{Z\}\. The two adaptive components are repeated until convergence or until the computational budget is exhausted\. A criterion based on the dominant variance contribution is proposed to switch between the two components at each iteration of the sequential algorithm\. This provides a principled framework for sequentially targeting the most prominent source of global predictive error across the input space𝒳\\mathcal\{X\}\.

The infilling of inducing inputs is discussed in[Section3\.1](https://arxiv.org/html/2606.31284#S3.SS1)\. The procedure for acquiring new observations of the data\-generating process is detailed in[Section3\.2](https://arxiv.org/html/2606.31284#S3.SS2)\. The switching criterion and the complete algorithm are presented in[Section3\.3](https://arxiv.org/html/2606.31284#S3.SS3)\. Computational considerations are discussed in[Section3\.4](https://arxiv.org/html/2606.31284#S3.SS4)\.

### 3\.1Inducing\-input infilling

In this section, we treat the training dataset𝒟\\mathcal\{D\}as fixed and address the infilling of the inducing\-input set𝒵\\mathcal\{Z\}\.

With fixed GP hyperparameters, adding new inducing inputs reduces the conditional\-prior variance\. Assuming that this addition does not significantly alter the remaining terms in[Eq\.27b](https://arxiv.org/html/2606.31284#S3.E27.2), global predictive accuracy can be enhanced by targeting the conditional\-prior variance contribution\. We thus propose to add inducing inputs to𝒵\\mathcal\{Z\}at locations that maximize the integrated reduction in conditional\-prior variance over the input space𝒳\\mathcal\{X\}\. This strategy is expected to yield performance comparable to, or exceeding that of, greedy strategies based on the maximum conditional\-prior variance, particularly in high\-dimensional settings\. See\[[26](https://arxiv.org/html/2606.31284#bib.bib26),[7](https://arxiv.org/html/2606.31284#bib.bib7),[2](https://arxiv.org/html/2606.31284#bib.bib2)\]for related work on optimal experimental design with standard GP regression\.

The conditional\-prior variance at an arbitrary input𝐱∈𝒳\\mathbf\{x\}\\in\\mathcal\{X\}, given the current inducing inputs𝐙\\mathbf\{Z\}and a candidate inducing input𝐳~∈𝒳\\tilde\{\\mathbf\{z\}\}\\in\\mathcal\{X\}, is expressed as

\(28\)Σz2​\(𝐱,𝐱;\[𝐙,𝐳~\]\)=K​\(𝐱,𝐱\)−\[K​\(𝐙,𝐱\)K​\(𝐳~,𝐱\)\]⊤​\[K​\(𝐙,𝐙\)K​\(𝐳~,𝐙\)⊤K​\(𝐳~,𝐙\)K​\(𝐳~,𝐳~\)\]−1​\[K​\(𝐙,𝐱\)K​\(𝐳~,𝐱\)\]\.\\Sigma\_\{z\}^\{2\}\(\\mathbf\{x\},\\mathbf\{x\}\\,;\[\\mathbf\{Z\},\\tilde\{\\mathbf\{z\}\}\]\)=K\(\\mathbf\{x\},\\mathbf\{x\}\)\-\\begin\{bmatrix\}K\(\\mathbf\{Z\},\\mathbf\{x\}\)\\\\ K\(\\tilde\{\\mathbf\{z\}\},\\mathbf\{x\}\)\\end\{bmatrix\}^\{\\top\}\\begin\{bmatrix\}K\(\\mathbf\{Z\},\\mathbf\{Z\}\)&K\(\\tilde\{\\mathbf\{z\}\},\\mathbf\{Z\}\)^\{\\top\}\\\\ K\(\\tilde\{\\mathbf\{z\}\},\\mathbf\{Z\}\)&K\(\\tilde\{\\mathbf\{z\}\},\\tilde\{\\mathbf\{z\}\}\)\\end\{bmatrix\}^\{\-1\}\\begin\{bmatrix\}K\(\\mathbf\{Z\},\\mathbf\{x\}\)\\\\ K\(\\tilde\{\\mathbf\{z\}\},\\mathbf\{x\}\)\\end\{bmatrix\}\.To avoid direct inversion of the block matrix, we apply the block matrix inversion identity based on the Schur complement, expressing the inverse in terms ofK​\(𝐙,𝐙\)−1K\(\\mathbf\{Z\},\\mathbf\{Z\}\)^\{\-1\}\. SinceK​\(𝐙,𝐙\)−1K\(\\mathbf\{Z\},\\mathbf\{Z\}\)^\{\-1\}is already available from the training of the quantile model,[Eq\.28](https://arxiv.org/html/2606.31284#S3.E28)becomes significantly cheaper to evaluate\. The integrated reduction in conditional\-prior variance induced by adding𝐳~\\tilde\{\\mathbf\{z\}\}is then

\(29a\)Δz​\(𝐳~;𝐙\)\\displaystyle\\Delta\_\{z\}\(\\tilde\{\\mathbf\{z\}\}\\,;\\mathbf\{Z\}\):=∫𝒳\{Σz2​\(𝐱,𝐱;𝐙\)−Σz2​\(𝐱,𝐱;\[𝐙​𝐳~\]\)\}​𝑑𝐱\\displaystyle:=\\int\_\{\\mathcal\{X\}\}\\Big\\\{\\Sigma\_\{z\}^\{2\}\(\\mathbf\{x\},\\mathbf\{x\}\\,;\\mathbf\{Z\}\)\-\\Sigma\_\{z\}^\{2\}\(\\mathbf\{x\},\\mathbf\{x\}\\,;\[\\mathbf\{Z\}\\ \\tilde\{\\mathbf\{z\}\}\]\)\\Big\\\}\\,d\\mathbf\{x\}\(29b\)=∫𝒳\{\(K​\(𝐳~,𝐱\)−K​\(𝐳~,𝐙\)​K​\(𝐙,𝐙\)−1​K​\(𝐙,𝐱\)\)2K​\(𝐳~,𝐳~\)−K​\(𝐳~,𝐙\)​K​\(𝐙,𝐙\)−1​K​\(𝐳~,𝐙\)⊤\}​𝑑𝐱\.\\displaystyle=\\int\_\{\\mathcal\{X\}\}\\Bigg\\\{\\frac\{\\left\(K\(\\tilde\{\\mathbf\{z\}\},\\mathbf\{x\}\)\-K\(\\tilde\{\\mathbf\{z\}\},\\mathbf\{Z\}\)K\(\\mathbf\{Z\},\\mathbf\{Z\}\)^\{\-1\}K\(\\mathbf\{Z\},\\mathbf\{x\}\)\\right\)^\{2\}\}\{K\(\\tilde\{\\mathbf\{z\}\},\\tilde\{\\mathbf\{z\}\}\)\-K\(\\tilde\{\\mathbf\{z\}\},\\mathbf\{Z\}\)K\(\\mathbf\{Z\},\\mathbf\{Z\}\)^\{\-1\}K\(\\tilde\{\\mathbf\{z\}\},\\mathbf\{Z\}\)^\{\\top\}\}\\Bigg\\\}\\,d\\mathbf\{x\}\.The integral over𝒳\\mathcal\{X\}can be estimated through Monte Carlo sampling\.

Throughout the infilling procedure, each new inducing input𝐳new∈𝒳\\mathbf\{z\}\_\{\\text\{new\}\}\\in\\mathcal\{X\}is selected as the maximizer of the integrated reduction in conditional\-prior variance \([29](https://arxiv.org/html/2606.31284#S3.E29)\):

\(30\)𝐳new=arg⁡max𝐳~∈𝒳⁡Δz​\(𝐳~;𝐙\)\.\\mathbf\{z\}\_\{\\text\{new\}\}=\\arg\\max\_\{\\tilde\{\\mathbf\{z\}\}\\in\\mathcal\{X\}\}\\Delta\_\{z\}\(\\tilde\{\\mathbf\{z\}\}\\,;\\mathbf\{Z\}\)\.The gradient of[Eq\.29b](https://arxiv.org/html/2606.31284#S3.E29.2)with respect to𝐳~\\tilde\{\\mathbf\{z\}\}is available in closed form\. Consequently, the optimization problem \([30](https://arxiv.org/html/2606.31284#S3.E30)\) can be solved with a deterministic, gradient\-based algorithm\.

### 3\.2Training\-data acquisition

We now address the problem of augmenting the training dataset𝒟\\mathcal\{D\}with additional observations of the stochastic processy​\(⋅\)y\(\\cdot\), while keeping the inducing\-input set𝒵\\mathcal\{Z\}fixed\.

Since the inducing inputs𝐙\\mathbf\{Z\}remain unchanged, enriching the training dataset does not significantly alter the conditional\-prior variance\. We thus propose to acquire observations at locations selected via rejection sampling, with acceptance probabilities governed by the posterior\-induced variance\. This adaptive strategy is expected to improve training\-data acquisition efficiency by focusing on regions where posterior uncertainty over the inducing variables induces high variability in the quantile predictions\.

The rejection sampling algorithm is based on an upper boundσref2\\sigma\_\{\\text\{ref\}\}^\{2\}for the posterior\-induced variance:

\(31\)σref2≥max𝐱∈𝒳⁡Σu2​\(𝐱,𝐱;𝐲,𝐗,𝐙\)\.\\sigma\_\{\\text\{ref\}\}^\{2\}\\geq\\max\_\{\\mathbf\{x\}\\in\\mathcal\{X\}\}\\Sigma\_\{u\}^\{2\}\(\\mathbf\{x\},\\mathbf\{x\}\\,;\\mathbf\{y\},\\mathbf\{X\},\\mathbf\{Z\}\)\.This upper bound serves as a reference value for the acceptance probabilities\. At each iteration of the rejection sampling, a candidate input location𝐱prop\\mathbf\{x\}\_\{\\text\{prop\}\}is drawn uniformly over the input space𝒳\\mathcal\{X\}and accepted if

\(32\)υ<faccept​\(Σu2​\(𝐱prop,𝐱prop;𝐲,𝐗,𝐙\)σref2\),\\upsilon<f\_\{\\text\{accept\}\}\\left\(\\frac\{\\Sigma\_\{u\}^\{2\}\(\\mathbf\{x\}\_\{\\text\{prop\}\},\\mathbf\{x\}\_\{\\text\{prop\}\}\\,;\\mathbf\{y\},\\mathbf\{X\},\\mathbf\{Z\}\)\}\{\\sigma\_\{\\text\{ref\}\}^\{2\}\}\\right\),whereυ\\upsilonis a random variable uniformly distributed in\(0,1\)\\left\(0,1\\right\), andfaccept:\[0,1\]→\[0,1\]f\_\{\\text\{accept\}\}:\[0,1\]\\to\[0,1\]is a monotone acceptance function\. In this work, we use the identity functionfaccept​\(x\)=xf\_\{\\text\{accept\}\}\(x\)=x\. Alternative acceptance functions can be used to adjust the selectivity of the sampling process\. Since the gradient of[Eq\.26](https://arxiv.org/html/2606.31284#S2.E26)with respect to𝐱\\mathbf\{x\}is available in closed form, we set the upper bound to the maximum ofΣu2​\(𝐱,𝐱;𝐲,𝐗,𝐙\)\\Sigma\_\{u\}^\{2\}\(\\mathbf\{x\},\\mathbf\{x\}\\,;\\mathbf\{y\},\\mathbf\{X\},\\mathbf\{Z\}\), found with a deterministic, gradient\-based algorithm\.

### 3\.3Sequential algorithm

[Section3\.1](https://arxiv.org/html/2606.31284#S3.SS1)addressed the inducing\-input infilling given a fixed training dataset𝒟\\mathcal\{D\}\.[Section3\.2](https://arxiv.org/html/2606.31284#S3.SS2)addressed the adaptive acquisition of new training data given a fixed inducing\-input set𝒵\\mathcal\{Z\}, and thus a fixed structural complexity of the posterior quantile model\. The present section combines them into a sequential algorithm that alternates between the two according to a switching criterion\.

#### 3\.3\.1Switching criterion

At each iteration, the sequential algorithm selects between adaptive training\-data acquisition and inducing\-input infilling based on the dominant variance contribution to the global error, as measured by the IMSE \([27](https://arxiv.org/html/2606.31284#S3.E27)\)\. A practical heuristic in this setting is to acquire new training data as long as the integrated posterior\-induced varianceΣu2​\(𝐱,𝐱;𝐲,𝐗,𝐙\)\\Sigma\_\{u\}^\{2\}\(\\mathbf\{x\},\\mathbf\{x\}\\,;\\mathbf\{y\},\\mathbf\{X\},\\mathbf\{Z\}\), evaluated over the input space𝒳\\mathcal\{X\}, dominates the integrated conditional\-prior varianceΣz2​\(𝐱,𝐱;𝐙\)\\Sigma\_\{z\}^\{2\}\(\\mathbf\{x\},\\mathbf\{x\}\\,;\\mathbf\{Z\}\), evaluated over the same domain\. Inducing\-input infilling is triggered once the reverse holds, formally when

\(33\)∫𝒳Σz2​\(𝐱,𝐱;𝐙\)​𝑑𝐱≥cratio​∫𝒳Σu2​\(𝐱,𝐱;𝐲,𝐗,𝐙\)​𝑑𝐱,\\int\_\{\\mathcal\{X\}\}\\Sigma\_\{z\}^\{2\}\(\\mathbf\{x\},\\mathbf\{x\}\\,;\\mathbf\{Z\}\)\\,d\\mathbf\{x\}\\geq c\_\{\\text\{ratio\}\}\\int\_\{\\mathcal\{X\}\}\\Sigma\_\{u\}^\{2\}\(\\mathbf\{x\},\\mathbf\{x\}\\,;\\mathbf\{y\},\\mathbf\{X\},\\mathbf\{Z\}\)\\,d\\mathbf\{x\},wherecratio\>0c\_\{\\text\{ratio\}\}\>0is a user\-defined constant\. This criterion implicitly determines the numberMMof inducing points required to adequately represent the posterior given the training dataset𝒟\\mathcal\{D\}\. The choicecratio=1c\_\{\\text\{ratio\}\}=1corresponds to balancing the two IMSE variance contributions\. This parameter may also be adjusted to reflect different computational costs associated with training\-data acquisition and inducing\-input infilling\. In the numerical experiments,cratio=1c\_\{\\text\{ratio\}\}=1yielded consistently satisfactory results and no further tuning was performed\. In practice, the integrals in[Eq\.33](https://arxiv.org/html/2606.31284#S3.E33)are estimated by Monte Carlo sampling\.

#### 3\.3\.2Pseudocode

The sequential algorithm is summarized in[Algorithm1](https://arxiv.org/html/2606.31284#alg1)\. Each update of the posterior quantile model involves solving the bilevel optimization problem \([21](https://arxiv.org/html/2606.31284#S2.E21)\)\. The rejection sampling step is embarrassingly parallel, since candidate locations are drawn independently\.

Algorithm 1Sequential sparse Gaussian process quantile regressionDefine initial training values

𝐲\(0\)\\mathbf\{y\}^\{\(0\)\}, initial training inputs

𝐗\(0\)\\mathbf\{X\}^\{\(0\)\}, initial inducing points

𝐙\(0\)\\mathbf\{Z\}^\{\(0\)\}\.

Train the posterior quantile model

qτ∣𝐲\(0\),𝐗\(0\),𝐙\(0\)q\_\{\\tau\}\\mid\\mathbf\{y\}^\{\(0\)\},\\mathbf\{X\}^\{\(0\)\},\\mathbf\{Z\}^\{\(0\)\}\.

Set the iteration counters:

S←0S\\leftarrow 0,

T←0T\\leftarrow 0\.

whilecomputational budget not exhausteddo

ifswitching criterion not satisfiedthen

Set the batch size

BB\.

Sample

BBtraining input locations

𝐗new∈𝒳B\\mathbf\{X\}\_\{\\text\{new\}\}\\in\\mathcal\{X\}^\{B\}adaptively\.

Acquire observations

𝐲new∈ℝB\\mathbf\{y\}\_\{\\text\{new\}\}\\in\\mathbb\{R\}^\{B\}of the stochastic process

y​\(⋅\)y\(\\cdot\)at

𝐗new\\mathbf\{X\}\_\{\\text\{new\}\}\.

Augment the training dataset:

𝐲\(S\+1\)←\[𝐲\(S\);𝐲new\]\\mathbf\{y\}^\{\(S\+1\)\}\\leftarrow\[\\mathbf\{y\}^\{\(S\)\};\\mathbf\{y\}\_\{\\text\{new\}\}\],

𝐗\(S\+1\)←\[𝐗\(S\),𝐗new\]\\mathbf\{X\}^\{\(S\+1\)\}\\leftarrow\[\\mathbf\{X\}^\{\(S\)\},\\mathbf\{X\}\_\{\\text\{new\}\}\]\.

Update the posterior quantile model to

qτ∣𝐲\(S\+1\),𝐗\(S\+1\),𝐙\(T\)q\_\{\\tau\}\\mid\\mathbf\{y\}^\{\(S\+1\)\},\\mathbf\{X\}^\{\(S\+1\)\},\\mathbf\{Z\}^\{\(T\)\}\.

Increment the iteration counter:

S←S\+1S\\leftarrow S\+1\.

else

Solve the infilling criterion for

𝐳new∈𝒳\\mathbf\{z\}\_\{\\text\{new\}\}\\in\\mathcal\{X\}\.

Augment the inducing\-input set:

𝐙\(T\+1\)←\[𝐙\(T\),𝐳new\]\\mathbf\{Z\}^\{\(T\+1\)\}\\leftarrow\[\\mathbf\{Z\}^\{\(T\)\},\\mathbf\{z\}\_\{\\text\{new\}\}\]\.

Update the posterior quantile model to

qτ∣𝐲\(S\),𝐗\(S\),𝐙\(T\+1\)q\_\{\\tau\}\\mid\\mathbf\{y\}^\{\(S\)\},\\mathbf\{X\}^\{\(S\)\},\\mathbf\{Z\}^\{\(T\+1\)\}\.

Increment the iteration counter:

T←T\+1T\\leftarrow T\+1\.

endif

endwhile

returnThe posterior quantile model

qτ∣𝐲\(S\),𝐗\(S\),𝐙\(T\)q\_\{\\tau\}\\mid\\mathbf\{y\}^\{\(S\)\},\\mathbf\{X\}^\{\(S\)\},\\mathbf\{Z\}^\{\(T\)\}

### 3\.4Computational considerations

Two distinct aspects govern the computational complexity of the proposed method: the cost of training the posterior quantile model for a fixed inducing\-input set and the cost of adapting the inducing representation itself\.

For a fixed inducing\-input set𝒵\\mathcal\{Z\}, the inference problem reduces to the optimization of the surrogate log\-posterior𝒥~N\\widetilde\{\\mathcal\{J\}\}\_\{N\}introduced in[Section2\.3](https://arxiv.org/html/2606.31284#S2.SS3)\. In contrast to sparse variational GP formulations, the dominant optimization variables consist of the inducing variables𝐮∈ℝM\\mathbf\{u\}\\in\\mathbb\{R\}^\{M\}and the model hyperparameters\. The posterior covariance is not optimized directly but is instead recovered from the inverse of the negative Hessian of the log\-posterior surrogate at the MAP estimate\. As a result, the number of dominant optimization variables associated with the inducing representation reduces from𝒪​\(M2\)\\mathcal\{O\}\(M^\{2\}\)\(a mean vector and a covariance matrix\) to𝒪​\(M\)\\mathcal\{O\}\(M\)\.

The availability of closed\-form expressions for the surrogate objective, its gradient, and its Hessian further enables the use of deterministic gradient\-based and second\-order optimization methods\. For fixed hyperparameters, evaluating the objective and its gradient scales essentially linearly with the numberNNof observations and the numberMMof inducing variables, while assembling the Hessian requires𝒪​\(N​M2\)\\mathcal\{O\}\(NM^\{2\}\)operations\. The subsequent factorization or inversion of the Hessian contributes an additional𝒪​\(M3\)\\mathcal\{O\}\(M^\{3\}\)cost\. Consequently, the dominant cost of a complete training step is approximately𝒪​\(N​M2\+M3\)\\mathcal\{O\}\(NM^\{2\}\+M^\{3\}\), up to constants depending on the optimization procedure and kernel evaluations\. Although the asymptotic cost per training iteration remains comparable to that of other sparse GP formulations, the reduced optimization space and the availability of exact derivatives typically result in faster convergence of the optimization procedure in practice\.

An additional distinction with respect to existing sparse GP approaches to quantile regression concerns the treatment of the inducing inputs\. Existing methodologies specify the numberMMof inducing points*a priori*and subsequently determine their locations using predefined allocation strategies, such as random sampling, quasi\-Monte Carlo \(QMC\) sequences, or clustering procedures includingkk\-means\.

In contrast, the present methodology treats the inducing\-input set as an adaptive component of the model\. The variance\-reduction criterion introduced in[Section3\.1](https://arxiv.org/html/2606.31284#S3.SS1)sequentially enriches the inducing representation by identifying locations that maximize the reduction of the integrated conditional\-prior variance\. Combined with the switching criterion \([33](https://arxiv.org/html/2606.31284#S3.E33)\), this procedure simultaneously determines both the number and the locations of the inducing points\. Model complexity is therefore increased only when justified by the remaining predictive uncertainty\. This adaptive construction seeks a parsimonious inducing representation that remains commensurate with the amount of information contained in the available training data, thereby avoiding the need to prescribe the number of inducing points beforehand\. In particular, the sequential strategy avoids the common situation in which an excessively large number of inducing points is selected*a priori*to guarantee predictive accuracy\. Instead, the inducing representation is refined only when the conditional\-prior variance becomes the dominant source of predictive uncertainty\.

## 4Numerical experiments

### 4\.1Benchmark functions

The performance and robustness of the proposed quantile regression method are assessed on two benchmark functions drawn from the literature\. These functions present a range of structural complexities, including nonlinearity, multimodality, and heteroskedasticity, which are characteristic of real\-world modeling challenges\.

- •Sabater\. The first benchmark is the function introduced by Sabater et al\.\[[25](https://arxiv.org/html/2606.31284#bib.bib25)\]to evaluate their Bayesian quantile optimization framework\. The stochastic process is defined as \(34\)f​\(𝐱,ω\)=∑d=1Df​\(xd,ω\),f\(\\mathbf\{x\},\\omega\)=\\sum\_\{d=1\}^\{D\}f\(\\mathrm\{x\}\_\{d\},\\omega\),where𝐱∈\[a,b\]D\\mathbf\{x\}\\in\[a,b\]^\{D\}witha=2a=2andb=8b=8,ω∈Ω\\omega\\in\\varOmegadenotes an outcome, and \(35\)f​\(x,ω\)=3−4​exp⁡\(−4​\(x−4\)2\)−5\.2​exp⁡\(−4​\(x−6\)2\)\+x−ab−a​ξ1​\(ω\)\+b−xb−a​ξ2​\(ω\)\.f\(\\mathrm\{x\},\\omega\)=3\-4\\exp\\left\(\-4\(\\mathrm\{x\}\-4\)^\{2\}\\right\)\-5\.2\\exp\\left\(\-4\(\\mathrm\{x\}\-6\)^\{2\}\\right\)\+\\frac\{\\mathrm\{x\}\-a\}\{b\-a\}\\,\\xi\_\{1\}\(\\omega\)\+\\frac\{b\-\\mathrm\{x\}\}\{b\-a\}\\,\\xi\_\{2\}\(\\omega\)\.The random variablesξ1:Ω→ℝ\\xi\_\{1\}:\\varOmega\\to\\mathbb\{R\}andξ2:Ω→ℝ\\xi\_\{2\}:\\varOmega\\to\\mathbb\{R\}are independent and distributed as \(36\)ξ1∼𝒰​\(0,10\),ξ2∼𝒩​\(1\.01,0\.71\)\.\\xi\_\{1\}\\sim\\mathcal\{U\}\\big\(0,10\\big\),\\quad\\xi\_\{2\}\\sim\\mathcal\{N\}\\big\(1\.01,0\.71\\big\)\.The stochastic process exhibits mild heteroskedasticity and asymmetric noise, resulting from the distinct distributions of the two independent random variables\. Consistent with the original setup, we evaluate performance at the0\.80\.8quantile level\. The0\.80\.8\-quantile function has2D2^\{D\}local minima\.
- •Michalewicz\. The second benchmark is a one\-dimensional function used by Torossian et al\.\[[31](https://arxiv.org/html/2606.31284#bib.bib31)\]to evaluate quantile regression methods\. The stochastic process, based on the Michalewicz function\[[6](https://arxiv.org/html/2606.31284#bib.bib6)\], is defined as \(37\)f​\(x,ω\)=−2​sin⁡\(x\)​sin30⁡\(x2π\)−0\.1​cos3⁡\(π​x10\)\|2−2​sin⁡\(x\)​sin30⁡\(x2π\)\|​ξ2​\(ω\),f\(\\mathrm\{x\},\\omega\)=\-2\\sin\(\\mathrm\{x\}\)\\sin^\{30\}\\left\(\\frac\{\\mathrm\{x\}^\{2\}\}\{\\pi\}\\right\)\-\\frac\{0\.1\\cos^\{3\}\\left\(\\frac\{\\pi\\mathrm\{x\}\}\{10\}\\right\)\}\{\\left\\lvert 2\-2\\sin\(\\mathrm\{x\}\)\\sin^\{30\}\\left\(\\frac\{\\mathrm\{x\}^\{2\}\}\{\\pi\}\\right\)\\right\\rvert\}\\,\\xi^\{2\}\(\\omega\),wherex∈\[0,4\]\\mathrm\{x\}\\in\[0,4\],ω∈Ω\\omega\\in\\varOmega, andξ:Ω→ℝ\\xi:\\varOmega\\to\\mathbb\{R\}is given by \(38\)ξ​\(ω\)=3​η​\(ω\)​1\{η​\(ω\)<0\}\+6​η​\(ω\)​1\{η​\(ω\)≥0\},withη∼𝒩​\(0,1\)\.\\xi\(\\omega\)=3\\eta\(\\omega\)\\,\\mathds\{1\}\_\{\\\{\\eta\(\\omega\)<0\\\}\}\+6\\eta\(\\omega\)\\,\\mathds\{1\}\_\{\\\{\\eta\(\\omega\)\\geq 0\\\}\},\\quad\\text\{with\}\\quad\\eta\\sim\\mathcal\{N\}\\big\(0,1\\big\)\.The stochastic process is characterized by pronounced heteroskedasticity and an asymmetric noise structure\. Performance is evaluated at the0\.90\.9quantile level, one of the two extreme levels considered in the original study\. The0\.90\.9\-quantile function exhibits sharply localized shape variations and is nearly flat elsewhere\.

### 4\.2Setup and reporting conventions

The training input locations𝐗\\mathbf\{X\}are systematically initialized uniformly over the input space\. The inducing inputs𝐙\\mathbf\{Z\}are initialized via a Halton sequence, a QMC sampling method\.

From[Section4\.4](https://arxiv.org/html/2606.31284#S4.SS4)onward, the experiments are repeated over6464independent training datasets\. Solid, dashed, and dotted lines indicate averages over the independent runs, while shaded regions represent the corresponding95%95\\%confidence intervals\. Moreover, the true conditionalτ\\tau\-quantiles of the synthetic benchmarks can be derived analytically\. We exploit them to normalize the reported results by

\(39\)ϵτ,true=∫𝒳\(qτ,true​\(𝐱\)−q¯τ,true\)2​𝑑𝐱,\\epsilon\_\{\{\\tau,\\text\{true\}\}\}=\\int\_\{\\mathcal\{X\}\}\\big\(q\_\{\\tau,\\text\{true\}\}\(\\mathbf\{x\}\)\-\\bar\{q\}\_\{\\tau,\\text\{true\}\}\\big\)^\{2\}\\,d\\mathbf\{x\},whereq¯τ,true\\bar\{q\}\_\{\\tau,\\text\{true\}\}denotes the population mean\.

### 4\.3Empirical validity of the Laplace approximation

This section examines the empirical validity of the Laplace approximation to the posterior distribution over the inducing variables\. We investigate the convergence of the posteriorp​\(𝐮∣𝐙,𝐲,𝐗\)p\(\\mathbf\{u\}\\mid\\mathbf\{Z\},\\mathbf\{y\},\\mathbf\{X\}\)to a Gaussian centered at the MAP estimate𝐮^\\hat\{\\mathbf\{u\}\}, with covariance given by the inverse of the negative log\-posterior Hessian𝐇N\\mathbf\{H\}\_\{N\}\. The assessment is carried out on the one\-dimensional variant of theSabaterfunction withM=10M=10inducing inputs\. Three training datasets of sizesN=100N=100,1,0001\{,\}000, and10,00010\{,\}000observations are considered, representing small, intermediate, and large training\-dataset regimes for this benchmark\.

Our analysis proceeds in two steps\. First, we examine how the posterior shape evolves with dataset size and whether it converges to the Laplace approximation\. We compare the ordered statistics of10510^\{5\}samples drawn from the true posteriorp​\(𝐮∣𝐙,𝐲,𝐗\)p\(\\mathbf\{u\}\\mid\\mathbf\{Z\},\\mathbf\{y\},\\mathbf\{X\}\)using MCMC against the same number drawn from its Laplace approximationp^​\(𝐮∣𝐙,𝐲,𝐗\)\\hat\{p\}\(\\mathbf\{u\}\\mid\\mathbf\{Z\},\\mathbf\{y\},\\mathbf\{X\}\), for each of the three dataset sizes\. The resulting Q–Q plots, one per componentm=1,…,Mm=1,\\ldots,Mof the multivariate distribution, are reported in[Fig\.1](https://arxiv.org/html/2606.31284#S4.F1)\.

![Refer to caption](https://arxiv.org/html/2606.31284v1/figures/validation-mcmc-laplace/sabater_qq_plot.png)Figure 1:Q–Q plots comparing the ordered statistics of samples drawn from the true posteriorp​\(𝐮∣𝐙,𝐲,𝐗\)p\(\\mathbf\{u\}\\mid\\mathbf\{Z\},\\mathbf\{y\},\\mathbf\{X\}\)through Markov chain Monte Carlo \(MCMC\), against those drawn from its Laplace approximation, for each componentm=1,…,10m=1,\\ldots,10of the inducing variables\. Results are shown for theSabater 1Dbenchmark function, using datasets withN=100N=100\(left\),1,0001\{,\}000\(middle\), and10,00010\{,\}000observations \(right\)\.As[Fig\.1](https://arxiv.org/html/2606.31284#S4.F1)shows, agreement between the two sets of ordered statistics improves consistently across all components as the dataset sizeNNincreases\. This confirms that the true posterior becomes increasingly well approximated by a Gaussian distribution\. For the largest dataset \(N=10,000N=10\{,\}000\), the two distributions appear nearly indistinguishable, while for the intermediate dataset \(N=1,000N=1\{,\}000\), the approximation is already reasonably accurate\. Noticeable discrepancies remain only for the smallest dataset \(N=100N=100\), which suggests that the Laplace approximation is less reliable in low\-data regimes\.

While the Q–Q plots capture marginal agreement in the inducing variables, they do not reflect potential discrepancies in the joint structure\. Thus, the second analysis compares the posterior\-induced varianceΣu2​\(𝐱,𝐱;𝐲,𝐗,𝐙\)\\Sigma\_\{u\}^\{2\}\(\\mathbf\{x\},\\mathbf\{x\}\\,;\\mathbf\{y\},\\mathbf\{X\},\\mathbf\{Z\}\)across the input space𝒳\\mathcal\{X\}, using the covariance matrix𝐂^\\hat\{\\mathbf\{C\}\}obtained either from the Laplace approximation or estimated empirically from MCMC samples\. The GP hyperparameters𝝍\\boldsymbol\{\\psi\}are fixed to the values obtained during training, that is, at the solution of the bilevel optimization problem \([21](https://arxiv.org/html/2606.31284#S2.E21)\)\. Predictions are evaluated on a uniform grid of10,00010\{,\}000points\. Results are shown in[Fig\.2](https://arxiv.org/html/2606.31284#S4.F2)\.

![Refer to caption](https://arxiv.org/html/2606.31284v1/figures/validation-mcmc-laplace/sabater_predictions.png)Figure 2:Comparison of posterior\-induced variance predictions with covariance obtained through Markov chain Monte Carlo \(MCMC\) versus using the Laplace approximation\. Results are shown for theSabater 1Dbenchmark function, using datasets withN=100N=100\(left\),1,0001\{,\}000\(middle\), and10,00010\{,\}000observations \(right\)\.[Figure2](https://arxiv.org/html/2606.31284#S4.F2)illustrates that the posterior\-induced variance predictions align increasingly well with those from MCMC asNNgrows, which is consistent with the Q–Q plots and further supports the asymptotic validity of the Laplace approximation\.

### 4\.4Inducing\-input infilling

With a sufficiently large and well\-chosen inducing\-input set𝒵\\mathcal\{Z\}, quantile predictions should approach the exact solution, provided that enough training data are available\. This section examines how our variance\-based infilling strategy \(IVR\) compares to predefined inducing\-input allocation strategies in this regime\.

QMC sampling via Halton sequences serves as the reference predefined inducing\-input allocation strategy\. Both infilling strategies are first evaluated on the one\-dimensional version of theSabaterfunction withN=10,000N=10\{,\}000observations\. The inducing inputs are initialized toM=10M=10to give the initial quantile model sufficient structural complexity to guide the adaptive infilling\.[Figure3\(a\)](https://arxiv.org/html/2606.31284#S4.F3.sf1)tracks the evolution of the conditional\-prior varianceΣz2​\(𝐱,𝐱;𝐙\)\\Sigma\_\{z\}^\{2\}\(\\mathbf\{x\},\\mathbf\{x\}\\,;\\mathbf\{Z\}\), integrated over the input space𝒳\\mathcal\{X\}, as each infilling strategy progresses\.[Figure3\(b\)](https://arxiv.org/html/2606.31284#S4.F3.sf2)illustrates how this is reflected in predictive accuracy, as measured by the IMSE\.

![Refer to caption](https://arxiv.org/html/2606.31284v1/x1.png)\(a\)Integrated conditional\-prior variance\.
![Refer to caption](https://arxiv.org/html/2606.31284v1/x2.png)\(b\)Integrated mean squared error\.

Figure 3:Evolution of the normalized \(a\) integrated conditional\-prior variance and \(b\) integrated mean squared error \(IMSE\) as functions of the numberMMof inducing points, for the adaptive inducing\-input infilling strategy \(IVR, black\) and quasi\-Monte Carlo sampling \(QMC, blue\)\. Results are shown for theSabater 1Dbenchmark function\.TheSabater 1Dfunction is one\-dimensional, smooth, and moderately multimodal\. As a result, only a few inducing points are needed to capture its structure\. As shown in[Fig\.3\(a\)](https://arxiv.org/html/2606.31284#S4.F3.sf1), the integrated conditional\-prior variance decays exponentially under adaptive infilling, outpacing that achieved by the Halton sequence\. This translates into a corresponding IMSE reduction, as seen in[Fig\.3\(b\)](https://arxiv.org/html/2606.31284#S4.F3.sf2), until the error plateaus nearM=17M=17inducing points\. On average over the6464independent runs, the adaptive infilling strategy outperforms the Halton sequence in terms of IMSE up to approximatelyM=31M=31, beyond which both strategies achieve comparable predictive accuracy\.

Notably, the integrated conditional\-prior variance continues to decay exponentially even after the IMSE has leveled off\. This illustrates that the adaptive infilling strategy remains robust once further reductions in conditional\-prior variance no longer improve predictions\. Such a decoupling suggests that the remaining predictive error is dominated by sources that are not addressed by the infilling criterion\.

[Figure4](https://arxiv.org/html/2606.31284#S4.F4)examines this transition by illustrating how the two IMSE variance components evolve under adaptive infilling\. The tracked quantities are the conditional\-prior varianceΣz2​\(𝐱,𝐱;𝐙\)\\Sigma\_\{z\}^\{2\}\(\\mathbf\{x\},\\mathbf\{x\}\\,;\\mathbf\{Z\}\)and the posterior\-induced varianceΣu2​\(𝐱,𝐱;𝐲,𝐗,𝐙\)\\Sigma\_\{u\}^\{2\}\(\\mathbf\{x\},\\mathbf\{x\}\\,;\\mathbf\{y\},\\mathbf\{X\},\\mathbf\{Z\}\), each integrated over the input space𝒳\\mathcal\{X\}\. Both curves intersect nearM=17M=17inducing points, coinciding with the onset of the IMSE plateau in[Fig\.3\(b\)](https://arxiv.org/html/2606.31284#S4.F3.sf2)\. Beyond this threshold, the variance induced by the posterior uncertainty of the inducing variables emerges as the dominant error contribution\. It would then be beneficial to stop the inducing\-input infilling and trigger data acquisition\.

![Refer to caption](https://arxiv.org/html/2606.31284v1/x3.png)Figure 4:Evolution of the normalized integrated conditional\-prior predictive variance \(solid line\) and integrated posterior\-induced variance \(dashed line\) as functions of the numberMMof inducing points\. Results are shown for theSabater 1Dbenchmark function\.The integrated variance terms analyzed here correspond precisely to the left\- and right\-hand sides of the switching criterion \([33](https://arxiv.org/html/2606.31284#S3.E33)\)\. As a result, the proposed switching criterion would halt the inducing\-input infilling within the regime of significant IMSE reduction and trigger data acquisition accordingly\. It proved effective in all our experiments\.

We now examine the inducing\-input infilling in higher dimension, reproducing the above experiments on the two\-dimensional variant of theSabaterfunction\. The quantile model is initialized withM=50M=50inducing points andN=10,000N=10\{,\}000observations\. The integrated conditional\-prior variance and the IMSE are reported in[Fig\.5](https://arxiv.org/html/2606.31284#S4.F5)\.

![Refer to caption](https://arxiv.org/html/2606.31284v1/x4.png)\(a\)Integrated conditional\-prior variance\.
![Refer to caption](https://arxiv.org/html/2606.31284v1/x5.png)\(b\)Integrated mean squared error\.

Figure 5:Evolution of the normalized \(a\) integrated conditional\-prior variance and \(b\) integrated mean squared error \(IMSE\) as functions of the numberMMof inducing points, for the adaptive inducing\-input infilling strategy \(IVR, black\) and quasi\-Monte Carlo sampling \(QMC, blue\)\. Results are shown for theSabater 2Dbenchmark function\.In line with the earlier findings, IVR consistently outperforms the predefined inducing\-input allocation strategy in terms of integrated conditional\-prior variance\. The exponential decay observed in the one\-dimensional case also holds in higher dimension\. In this example, the IMSE achieved through IVR does not plateau but instead reaches a minimum before rising again\. The minimum is reached nearM=140M=140inducing points\. Compared with the one\-dimensional case, this indicates that more inducing points are needed to capture the increased structural complexity of higher\-dimensional settings\. The IMSE deterioration occurs because the available training data become insufficient to reliably estimate the inducing variables, causing the posterior\-induced variance contribution to increase\. As established previously, the switching criterion \([33](https://arxiv.org/html/2606.31284#S3.E33)\) would guard against such a deterioration\.

### 4\.5Predictive accuracy with different training dataset sizes

The previous experiments showed that, provided sufficient training data, our adaptive inducing\-input infilling strategy outperforms predefined allocation strategies\. It remains well\-behaved throughout the infilling, even when the targeted variance contribution in the IMSE becomes dominated by other contributions\. This section investigates how the global predictive accuracy of the quantile model varies with the numberNNof observations in the training dataset𝒟\\mathcal\{D\}\. We study how the optimal number of inducing points evolves with increasing data availability, and how the model behaves under limited data regimes\.

Three training datasets of sizeN=300N=300,3,0003\{,\}000, and30,00030\{,\}000are considered for the one\-dimensionalSabaterbenchmark\. We compare the evolution of the IMSE throughout the adaptive inducing\-input infilling for each dataset size, starting fromM=10M=10inducing points\. Results are reported in[Fig\.6](https://arxiv.org/html/2606.31284#S4.F6)\.

![Refer to caption](https://arxiv.org/html/2606.31284v1/x6.png)Figure 6:Evolution of the integrated mean squared error \(IMSE\) as a function of the numberMMof inducing points, for training dataset sizesN=300N=300\(dotted line\),3,0003\{,\}000\(dashed line\), and30,00030\{,\}000\(solid line\)\. Results are shown for theSabater 1Dbenchmark function\.As shown in[Fig\.6](https://arxiv.org/html/2606.31284#S4.F6), the number of inducing points required to reach the minimum IMSE increases with the number of observations\. This indicates that larger datasets can support more complex inducing representations\. Moreover, the minimum IMSE decreases asNNgrows, confirming that the best achievable predictive accuracy improves with more data\. In the asymptotic regimeN⟶∞N\\longrightarrow\\infty, the complexity of the trueτ\\tau\-quantile function can be fully captured, and predictions should converge to the exact values\. These observations suggest that the proposed strategy is compatible with asymptotic convergence toward the true quantile function as both the training set size and the inducing\-point complexity increase\.

With limited training data \(N=300N=300andN=3,000N=3\{,\}000\), the IMSE does not plateau but instead reaches a minimum before rising again, for the same reason as in the two\-dimensional example of[Section4\.4](https://arxiv.org/html/2606.31284#S4.SS4)\. This further illustrates the relevance of the switching criterion \([33](https://arxiv.org/html/2606.31284#S3.E33)\)\.

### 4\.6Sequential enrichment

[Section4\.4](https://arxiv.org/html/2606.31284#S4.SS4)established that infilling the inducing\-input set by maximizing the integrated reduction in conditional\-prior variance outperforms predefined allocation strategies, provided sufficient training data\.[Section4\.5](https://arxiv.org/html/2606.31284#S4.SS5)showed that global predictive accuracy improves with training dataset size, and that the optimal number of inducing points grows accordingly\. The variance\-based switching criterion \([33](https://arxiv.org/html/2606.31284#S3.E33)\) is designed to identify the dominant variance contribution in the IMSE, and to trigger enrichment of either the inducing\-input set𝒵\\mathcal\{Z\}or the training dataset𝒟\\mathcal\{D\}accordingly\. This section compares the convergence of the sequential algorithm, which comprises inducing\-input infilling and training\-data acquisition, as a function of training dataset size, for two acquisition strategies: uniform sampling and the rejection sampling scheme introduced in[Section3\.2](https://arxiv.org/html/2606.31284#S3.SS2)\.

The comparison is carried out on both theSabater 2DandMichalewiczfunctions\. For theSabater 2Dbenchmark, the initial number of training data and inducing points are respectivelyN=500N=500andM=50M=50\. TheMichalewiczbenchmark is initialized withN=300N=300observations andM=10M=10inducing points\. In both cases, subsequent inducing inputs are added through IVR\. Training data are acquired in batches of size equal to10%10\\%of the current dataset\. Results are reported in[Fig\.7](https://arxiv.org/html/2606.31284#S4.F7)\.

For theSabater 2Dbenchmark \([Fig\.7\(a\)](https://arxiv.org/html/2606.31284#S4.F7.sf1)\), rejection sampling performs on par with uniform sampling\. This suggests that, for this function, the benefits of adaptivity are concentrated in the inducing\-input allocation rather than in the data\-acquisition strategy\. In contrast, for theMichalewiczbenchmark \([Fig\.7\(b\)](https://arxiv.org/html/2606.31284#S4.F7.sf2)\), rejection sampling yields substantially lower IMSE values than uniform sampling up to approximately2,0002\{,\}000observations, showing that adaptive data acquisition can provide significant gains over predefined strategies, even for one\-dimensional functions\. TheMichalewiczbenchmark exhibits localized shape variations that rejection sampling helps to discover early\. This early discovery then guides the inference toward predictions that significantly outperform those obtained with uniform sampling, until the latter has acquired enough data for the model to recover the predictive accuracy of the rejection\-based model\. We therefore conclude that adaptive data acquisition is particularly beneficial in data\-limited settings for functions with localized shape variations\. Overall, in both cases, the IMSE decreases consistently as the algorithm progresses, supporting the argument of[Section4\.5](https://arxiv.org/html/2606.31284#S4.SS5)that the sequential adaptive strategy converges asymptotically to the trueτ\\tau\-quantile function\.

![Refer to caption](https://arxiv.org/html/2606.31284v1/x7.png)\(a\)Sabater 2D
![Refer to caption](https://arxiv.org/html/2606.31284v1/x8.png)\(b\)Michalewicz

Figure 7:Evolution of the integrated mean squared error \(IMSE\) as a function of the numberNNof observations, for the sequential algorithm with data acquired through rejection sampling \(black\) and uniform sampling \(red\)\. Results are shown for \(a\) theSabater 2Dand \(b\) theMichalewiczbenchmark functions\.

## 5Conclusion

This work introduced a sequential sparse Gaussian process framework for Bayesian quantile regression\. The proposed approach combines a sparse representation of the latent quantile function through inducing variables with a Laplace approximation to their posterior distribution\. To avoid the repeated sampling that would otherwise be required to evaluate the exact posterior objective, inference is recast as the optimization of a tractable surrogate log\-posterior, yielding a sample\-free Bayesian inference procedure together with closed\-form expressions for the predictive mean, covariance, and their derivatives\. The resulting formulation leads to a reduced optimization space compared with sparse variational approaches and enables the use of efficient deterministic, gradient\-based optimization algorithms\.

A central outcome of this work is the decomposition of the predictive uncertainty into two complementary contributions: a conditional\-prior variance associated with the inducing representation and a posterior\-induced variance associated with uncertainty in the inducing variables\. This decomposition provides a principled framework for sequentially improving the predictive quantile model\. The conditional\-prior variance motivates an adaptive inducing\-input infilling strategy that progressively enriches the sparse representation of the quantile function\. The posterior\-induced variance, in turn, drives an adaptive data\-acquisition procedure that focuses in regions where uncertainty in the inducing variables has the largest impact on the quantile predictions\. These two enrichment mechanisms are unified through a variance\-based switching criterion that balances model\-complexity growth and data acquisition according to the dominant source of predictive uncertainty\.

The numerical experiments demonstrate several important properties of the proposed methodology\. First, the Laplace approximation provides accurate estimates of the posterior\-induced covariance when compared with reference Markov chain Monte Carlo computations\. Second, the variance\-based inducing\-input infilling consistently outperforms predefined allocation strategies in terms of conditional\-prior variance reduction and predictive accuracy\. Third, the evolution of the two variance contributions confirms the rationale of the proposed switching criterion, which identifies the regime where further enrichment of the inducing representation yields diminishing returns and additional data acquisition becomes preferable\. Finally, the complete sequential strategy improves predictive accuracy while adaptively controlling the complexity of the sparse representation\.

The present work also highlights several directions for future research\. From a theoretical perspective, a more detailed analysis of the surrogate objective and its impact on the location of the maximum*a posteriori*estimate would further strengthen the foundations of the approach\. Similarly, establishing convergence guarantees for the sequential enrichment procedure remains an open problem\. From a methodological perspective, the framework could be extended to heteroskedastic quantile models and simultaneous estimation of multiple quantile levels\. These developments constitute promising directions for further advancing Bayesian quantile regression under limited computational and data\-acquisition budgets\.

## Appendix AControl of the lower\-bound surrogate error

This appendix provides a local justification of the surrogate log\-joint density introduced in \([17](https://arxiv.org/html/2606.31284#S2.E17)\)\. We quantify the approximation error in terms of the conditional\-prior varianceΣz2​\(𝐱,𝐱;𝐙\)\\Sigma\_\{z\}^\{2\}\(\\mathbf\{x\},\\mathbf\{x\}\\,;\\mathbf\{Z\}\), which is the variance component directly targeted by the inducing\-input infilling procedure\. The following notation is introduced to simplify the expression of the GP prediction at the observations:

μn:=μ​\(𝐱n;𝐮,𝐙,𝝍\),δn∼𝒩​\(0,σn2\),σn:=Σz2​\(𝐱n,𝐱n;𝐙,𝝍\)\.\\mu\_\{n\}:=\\mu\(\\mathbf\{x\}\_\{n\}\\,;\\mathbf\{u\},\\mathbf\{Z\},\\boldsymbol\{\\psi\}\),\\qquad\\delta\_\{n\}\\sim\\mathcal\{N\}\\big\(0,\\sigma\_\{n\}^\{2\}\\big\),\\qquad\\sigma\_\{n\}:=\\sqrt\{\\Sigma\_\{z\}^\{2\}\(\\mathbf\{x\}\_\{n\},\\mathbf\{x\}\_\{n\}\\,;\\mathbf\{Z\},\\boldsymbol\{\\psi\}\)\}\.
For the asymmetric Laplace likelihood in[Eq\.4](https://arxiv.org/html/2606.31284#S2.E4), the log\-likelihood of a single training point\(𝐱n,yn\)\(\\mathbf\{x\}\_\{n\},\\mathrm\{y\}\_\{n\}\)is

ℓn​\(δn\)=log⁡\(τ​\(1−τ\)α\)−1α​ρτ​\(yn−μn−δn\)\.\\ell\_\{n\}\(\\delta\_\{n\}\)=\\log\\left\(\\frac\{\\tau\(1\-\\tau\)\}\{\\alpha\}\\right\)\-\\frac\{1\}\{\\alpha\}\\rho\_\{\\tau\}\(\\mathrm\{y\}\_\{n\}\-\\mu\_\{n\}\-\\delta\_\{n\}\)\.The exact and surrogate observation\-wise contributions are

jn:=log⁡\(𝔼​\[exp⁡\{ℓn​\(δn\)\}\]\),j~n:=𝔼​\[ℓn​\(δn\)\]\.j\_\{n\}:=\\log\\left\(\\mathbb\{E\}\\big\[\\exp\\\{\\ell\_\{n\}\(\\delta\_\{n\}\)\\\}\\big\]\\right\),\\qquad\\widetilde\{j\}\_\{n\}:=\\mathbb\{E\}\\big\[\\ell\_\{n\}\(\\delta\_\{n\}\)\\big\]\.We setℓ¯n:=ℓn−𝔼​\[ℓn\]\\bar\{\\ell\}\_\{n\}:=\\ell\_\{n\}\-\\mathbb\{E\}\\big\[\\ell\_\{n\}\\big\], the centered log\-likelihood of the observation, to recast the difference in the observation\-wise contributions as

jn−j~n=log⁡\(𝔼​\[exp⁡\(ℓ¯n\)\]\),j\_\{n\}\-\\widetilde\{j\}\_\{n\}=\\log\\left\(\\mathbb\{E\}\\big\[\\exp\(\\bar\{\\ell\}\_\{n\}\)\\big\]\\right\),which is nonnegative by Jensen’s inequality\. It also shows that this difference equals the cumulant generating functionKn​\(t\)K\_\{n\}\(t\)ofℓ¯n\\bar\{\\ell\}\_\{n\}evaluated att=1t=1\. Assuming thatKn​\(t\)K\_\{n\}\(t\)remains finite in a neighborhood of the origin, the cumulant generating function has the expansion

Kn​\(t\):=log⁡\(𝔼​\[exp⁡\(t​ℓ¯n\)\]\)=∑r≥2κr,nr\!​tr,K\_\{n\}\(t\):=\\log\\left\(\\mathbb\{E\}\\big\[\\exp\(t\\bar\{\\ell\}\_\{n\}\)\\big\]\\right\)=\\sum\_\{r\\geq 2\}\\frac\{\\kappa\_\{r,n\}\}\{r\!\}t^\{r\},whereκr,n\\kappa\_\{r,n\}is therr\-th cumulant ofℓn​\(δn\)\\ell\_\{n\}\(\\delta\_\{n\}\)\. Evaluating this expansion att=1t=1yields the formal expansion of the difference between the contributions

jn−j~n=12​𝕍​\[ℓn​\(δn\)\]\+16​κ3,n\+124​κ4,n\+⋯\.j\_\{n\}\-\\widetilde\{j\}\_\{n\}=\\frac\{1\}\{2\}\\mathbb\{V\}\\big\[\\ell\_\{n\}\(\\delta\_\{n\}\)\\big\]\+\\frac\{1\}\{6\}\\kappa\_\{3,n\}\+\\frac\{1\}\{24\}\\kappa\_\{4,n\}\+\\cdots\.Consequently, when the conditional\-prior varianceσn2\\sigma\_\{n\}^\{2\}is small, the leading term in the gap between the objectives is

jn−j~n=12​𝕍​\[ℓn​\(δn\)\]\+𝒪​\(𝔼​\[\|ℓ¯n\|3\]\)\.j\_\{n\}\-\\widetilde\{j\}\_\{n\}=\\frac\{1\}\{2\}\\mathbb\{V\}\\big\[\\ell\_\{n\}\(\\delta\_\{n\}\)\\big\]\+\\mathcal\{O\}\\left\(\\mathbb\{E\}\\big\[\\lvert\\bar\{\\ell\}\_\{n\}\\rvert^\{3\}\\big\]\\right\)\.This identity makes explicit that the surrogate is accurate whenever the log\-likelihood fluctuates weakly under the conditional GP prior\. Moreover, since the check function \([3](https://arxiv.org/html/2606.31284#S2.E3)\) is globally Lipschitz continuous with constantLτ:=max⁡\(τ,1−τ\)L\_\{\\tau\}:=\\max\(\\tau,1\-\\tau\), it follows that

\|ℓn​\(δ\)−ℓn​\(δ′\)\|≤Lτα​\|δ−δ′\|\.\\lvert\\ell\_\{n\}\(\\delta\)\-\\ell\_\{n\}\(\\delta^\{\\prime\}\)\\rvert\\leq\\frac\{L\_\{\\tau\}\}\{\\alpha\}\\lvert\\delta\-\\delta^\{\\prime\}\\rvert\.In particular, sinceδn∼𝒩​\(0,σn2\)\\delta\_\{n\}\\sim\\mathcal\{N\}\\big\(0,\\sigma\_\{n\}^\{2\}\\big\), the Gaussian Poincaré inequality gives

𝕍​\[ℓn​\(δn\)\]≤Lτ2α2​σn2\.\\mathbb\{V\}\\big\[\\ell\_\{n\}\(\\delta\_\{n\}\)\\big\]\\leq\\frac\{L\_\{\\tau\}^\{2\}\}\{\\alpha^\{2\}\}\\sigma\_\{n\}^\{2\}\.The third absolute central moment can be bounded analogously, sinceℓn\\ell\_\{n\}is a Lipschitz transform of a Gaussian random variable\. Thus, for smallσn\\sigma\_\{n\},

jn−j~n≤Lτ22​α2​σn2\+𝒪​\(σn3\)\.j\_\{n\}\-\\widetilde\{j\}\_\{n\}\\leq\\frac\{L\_\{\\tau\}^\{2\}\}\{2\\alpha^\{2\}\}\\sigma\_\{n\}^\{2\}\+\\mathcal\{O\}\\left\(\\sigma\_\{n\}^\{3\}\\right\)\.Summing over the observations gives the difference between the exact and surrogate conditional log\-likelihoods\. ReintroducingΣz2\\Sigma\_\{z\}^\{2\}into the notation and dividing byNNthen yields

0≤1N​\(𝒥N​\(𝐮\)−𝒥~N​\(𝐮\)\)≤Lτ22​α2​1N​∑n=1NΣz2​\(𝐱n,𝐱n;𝐙\)\+𝒪​\(1N​∑n=1N\(Σz2​\(𝐱n,𝐱n;𝐙\)\)3/2\)\.0\\leq\\frac\{1\}\{N\}\\left\(\{\\mathcal\{J\}\}\_\{N\}\(\\mathbf\{u\}\)\-\\widetilde\{\\mathcal\{J\}\}\_\{N\}\(\\mathbf\{u\}\)\\right\)\\leq\\frac\{L\_\{\\tau\}^\{2\}\}\{2\\alpha^\{2\}\}\\frac\{1\}\{N\}\\sum\_\{n=1\}^\{N\}\\Sigma\_\{z\}^\{2\}\(\\mathbf\{x\}\_\{n\},\\mathbf\{x\}\_\{n\}\\,;\\mathbf\{Z\}\)\+\\mathcal\{O\}\\left\(\\frac\{1\}\{N\}\\sum\_\{n=1\}^\{N\}\(\\Sigma\_\{z\}^\{2\}\(\\mathbf\{x\}\_\{n\},\\mathbf\{x\}\_\{n\}\\,;\\mathbf\{Z\}\)\)^\{3/2\}\\right\)\.Therefore, if the cumulant expansion is valid, the gap between the exact conditional log\-likelihood and the surrogate vanishes whenever the average conditional\-prior variance at the training inputs vanishes\. In particular, if the sequence of inducing sets𝒵\\mathcal\{Z\}is such that1N​∑n=1NΣz2​\(𝐱n,𝐱n;𝐙N\)⟶0,\\frac\{1\}\{N\}\\sum\_\{n=1\}^\{N\}\\Sigma\_\{z\}^\{2\}\(\\mathbf\{x\}\_\{n\},\\mathbf\{x\}\_\{n\}\\,;\\mathbf\{Z\}\_\{N\}\)\\longrightarrow 0,then1N​\(𝒥N​\(𝐮\)−𝒥~N​\(𝐮\)\)⟶0\.\\frac\{1\}\{N\}\\left\(\{\\mathcal\{J\}\}\_\{N\}\(\\mathbf\{u\}\)\-\\widetilde\{\\mathcal\{J\}\}\_\{N\}\(\\mathbf\{u\}\)\\right\)\\longrightarrow 0\.This result does not imply that increasingNNalone improves the surrogate\. If the inducing set𝒵\\mathcal\{Z\}is fixed, the conditional\-prior varianceΣz2​\(⋅,⋅;𝐙\)\\Sigma\_\{z\}^\{2\}\(\\cdot,\\cdot\\,;\\mathbf\{Z\}\)is fixed, and the total surrogate gap may grow linearly withNN\. The relevant asymptotic regime is therefore one in which the data size increases while the inducing set is enriched, so that the average conditional\-prior variance at the training inputs decreases to zero\. This is precisely the regime promoted by the inducing\-input infilling strategy\.

## Appendix BApproximate log\-joint density and its Hessian

Here, we provide the expressions of the approximate log\-joint density𝒥~N\\widetilde\{\\mathcal\{J\}\}\_\{N\}and its Hessian\. Using the notation introduced in[AppendixA](https://arxiv.org/html/2606.31284#A1), the approximate log\-joint density𝒥~N\\widetilde\{\\mathcal\{J\}\}\_\{N\}is given by

\(40\)𝒥~N​\(𝐮;𝝍,α\)\\displaystyle\\widetilde\{\\mathcal\{J\}\}\_\{N\}\(\\mathbf\{u\}\\,;\\boldsymbol\{\\psi\},\\alpha\)=𝔼π​\(qτ∣𝐮,𝐙,𝝍\)​\[log⁡p​\(𝐲∣𝐗,qτ,α\)\]\+log⁡π​\(𝐮∣𝐙,𝝍\)\\displaystyle=\\mathbb\{E\}\_\{\\pi\(q\_\{\\tau\}\\mid\\mathbf\{u\},\\mathbf\{Z\},\\boldsymbol\{\\psi\}\)\}\\big\[\\log p\(\\mathbf\{y\}\\mid\\mathbf\{X\},q\_\{\\tau\},\\alpha\)\\big\]\+\\log\\pi\(\\mathbf\{u\}\\mid\\mathbf\{Z\},\\boldsymbol\{\\psi\}\)=−1α​∑n=1N\{\(yn−μn\)​\[Φ​\(yn−μnσn\)\+τ−1\]\+σn​φ​\(yn−μnσn\)\}\\displaystyle=\-\\frac\{1\}\{\\alpha\}\\sum\_\{n=1\}^\{N\}\\Bigg\\\{\\left\(\\mathrm\{y\}\_\{n\}\-\\mu\_\{n\}\\right\)\\left\[\\Phi\\left\(\\frac\{\\mathrm\{y\}\_\{n\}\-\\mu\_\{n\}\}\{\\sigma\_\{n\}\}\\right\)\+\\tau\-1\\right\]\+\\sigma\_\{n\}\\,\\varphi\\left\(\\frac\{\\mathrm\{y\}\_\{n\}\-\\mu\_\{n\}\}\{\\sigma\_\{n\}\}\\right\)\\Bigg\\\}−N​log⁡α−12​𝐮⊤​K​\(𝐙,𝐙;𝝍\)​𝐮−12​log⁡\|K​\(𝐙,𝐙;𝝍\)\|\+Cst,\\displaystyle\\hskip 11\.9501pt\-N\\log\\alpha\-\\frac\{1\}\{2\}\\mathbf\{u\}^\{\\top\}K\(\\mathbf\{Z\},\\mathbf\{Z\}\\,;\\boldsymbol\{\\psi\}\)\\mathbf\{u\}\-\\frac\{1\}\{2\}\\log\\,\\left\\lvert K\(\\mathbf\{Z\},\\mathbf\{Z\}\\,;\\boldsymbol\{\\psi\}\)\\right\\rvert\+\\mathrm\{C^\{st\}\},whereΦ\\Phiandφ\\varphidenote the cumulative and probability density functions of the standard normal distribution, respectively\. The constant termCst\\mathrm\{C^\{st\}\}is independent of𝐮\\mathbf\{u\},𝝍\\boldsymbol\{\\psi\}, andα\\alpha\.

The Hessian of the approximate log\-joint density𝒥~N\\widetilde\{\\mathcal\{J\}\}\_\{N\}with respect to𝐮\\mathbf\{u\}is given by

\(41\)∇𝐮2𝒥~N​\(𝐮;𝝍,α\)\\displaystyle\\nabla\_\{\\mathbf\{u\}\}^\{2\}\\,\\widetilde\{\\mathcal\{J\}\}\_\{N\}\(\\mathbf\{u\}\\,;\\boldsymbol\{\\psi\},\\alpha\)=∇𝐮2𝔼π​\(qτ∣𝐮,𝐙,𝝍\)​\[log⁡p​\(𝐲∣𝐗,qτ,α\)\]\+∇𝐮2log⁡π​\(𝐮∣𝐙,𝝍\)\\displaystyle=\\nabla\_\{\\mathbf\{u\}\}^\{2\}\\,\\mathbb\{E\}\_\{\\pi\(q\_\{\\tau\}\\mid\\mathbf\{u\},\\mathbf\{Z\},\\boldsymbol\{\\psi\}\)\}\\big\[\\log p\(\\mathbf\{y\}\\mid\\mathbf\{X\},q\_\{\\tau\},\\alpha\)\\big\]\+\\nabla\_\{\\mathbf\{u\}\}^\{2\}\\log\\pi\(\\mathbf\{u\}\\mid\\mathbf\{Z\},\\boldsymbol\{\\psi\}\)=−1α∑n=1N\{K\(𝐱n,𝐙;𝝍\)K\(𝐙,𝐙;𝝍\)−1K\(𝐙,𝐙;𝝍\)−1K\(𝐱n,𝐙;𝝍\)⊤\\displaystyle=\-\\frac\{1\}\{\\alpha\}\\sum\_\{n=1\}^\{N\}\\Bigg\\\{K\(\\mathbf\{x\}\_\{n\},\\mathbf\{Z\}\\,;\\boldsymbol\{\\psi\}\)K\(\\mathbf\{Z\},\\mathbf\{Z\}\\,;\\boldsymbol\{\\psi\}\)^\{\-1\}K\(\\mathbf\{Z\},\\mathbf\{Z\}\\,;\\boldsymbol\{\\psi\}\)^\{\-1\}K\(\\mathbf\{x\}\_\{n\},\\mathbf\{Z\}\\,;\\boldsymbol\{\\psi\}\)^\{\\top\}×\[1σnφ\(yn−μnσn\)\]\}−K\(𝐙,𝐙;𝝍\)−1\.\\displaystyle\\hskip 61\.17325pt\\times\\left\[\\frac\{1\}\{\\sigma\_\{n\}\}\\,\\varphi\\left\(\\frac\{\\mathrm\{y\}\_\{n\}\-\\mu\_\{n\}\}\{\\sigma\_\{n\}\}\\right\)\\right\]\\Bigg\\\}\-K\(\\mathbf\{Z\},\\mathbf\{Z\}\\,;\\boldsymbol\{\\psi\}\)^\{\-1\}\.

## Acknowledgments

The work of the first author is supported by the LabCom MATritime funded by the Agence Nationale de la Recherche \(Grant No\. ANR\-22\-LCV2\-0010\)\. Numerical experiments presented in this paper were carried out using the PlaFRIM experimental testbed, supported by Inria, CNRS \(LaBRI and IMB\), Université de Bordeaux, Bordeaux INP, and Conseil Régional d’Aquitaine \(see https://www\.plafrim\.fr\)\.

## References

- \[1\]S\. Abeywardana and F\. Ramos,Variational inference for nonparametric Bayesian quantile regression, in Proceedings of the Twenty\-Ninth AAAI Conference on Artificial Intelligence, vol\. 29, Austin, Texas, USA, 2015,[https://doi\.org/10\.1609/aaai\.v29i1\.9465](https://doi.org/10.1609/aaai.v29i1.9465)\.
- \[2\]J\. Bect, D\. Ginsbourger, L\. Li, V\. Picheny, and E\. Vazquez,Sequential design of computer experiments for the estimation of a probability of failure, Statistics and Computing, 22 \(2012\), pp\. 773–793,[https://doi\.org/10\.1007/s11222\-011\-9241\-4](https://doi.org/10.1007/s11222-011-9241-4),[https://arxiv\.org/abs/1009\.5177](https://arxiv.org/abs/1009.5177)\.
- \[3\]A\. Boukouvalas, R\. Barillec, and D\. Cornford,Direct Gaussian process quantile regression using expectation propagation, in Proceedings of the 29th International Coference on Machine Learning, ICML’12, Madison, Wisconsin, USA, 2012, Omnipress, pp\. 939–946\.
- \[4\]D\. R\. Burt, C\. E\. Rasmussen, and M\. van der Wilk,Convergence of sparse variational inference in Gaussian processes regression, Journal of Machine Learning Research, 21 \(2020\), pp\. 1–63\.
- \[5\]A\. J\. Cannon,Quantile regression neural networks: implementation in R and application to precipitation downscaling, Computers & Geosciences, 37 \(2011\), pp\. 1277–1284,[https://doi\.org/10\.1016/j\.cageo\.2010\.07\.005](https://doi.org/10.1016/j.cageo.2010.07.005)\.
- \[6\]L\. C\. W\. Dixon and G\. P\. Szego,The global optimization problem: an introduction, Towards Global Optimiation 2, \(1978\), pp\. 1–15\.
- \[7\]R\. B\. Gramacy and H\. K\. H\. Lee,Adaptive design and analysis of supercomputer experiments, Technometrics, 51 \(2009\), pp\. 130–145,[https://doi\.org/10\.48550/arXiv\.0805\.4359](https://doi.org/10.48550/arXiv.0805.4359),[https://arxiv\.org/abs/0805\.4359](https://arxiv.org/abs/0805.4359)\.
- \[8\]J\. Hensman, N\. Fusi, and N\. D\. Lawrence,Gaussian processes for big data, 2013,[https://doi\.org/10\.48550/arXiv\.1309\.6835](https://doi.org/10.48550/arXiv.1309.6835),[https://arxiv\.org/abs/1309\.6835](https://arxiv.org/abs/1309.6835)\.
- \[9\]J\. Hensman, A\. G\. Matthews, M\. Filippone, and Z\. Ghahramani,MCMC for variationally sparse Gaussian processes, in Advances in Neural Information Processing Systems, vol\. 28, Curran Associates, Inc\., 2015\.
- \[10\]S\. R\. Jantre, S\. Bhattacharya, and T\. Maiti,Quantile regression neural networks: a Bayesian approach, Journal of Statistical Theory and Practice, 15 \(2021\), p\. 68,[https://doi\.org/10\.1007/s42519\-021\-00189\-w](https://doi.org/10.1007/s42519-021-00189-w),[https://arxiv\.org/abs/2009\.13591](https://arxiv.org/abs/2009.13591)\.
- \[11\]R\. Koenker and G\. Bassett,Regression quantiles, Econometrica, 46 \(1978\), pp\. 33–50,[https://doi\.org/10\.2307/1913643](https://doi.org/10.2307/1913643),[https://arxiv\.org/abs/1913643](https://arxiv.org/abs/1913643)\.
- \[12\]R\. Koenker, P\. Ng, and S\. Portnoy,Quantile smoothing splines, Biometrika, 81 \(1994\), pp\. 673–680,[https://doi\.org/10\.1093/biomet/81\.4\.673](https://doi.org/10.1093/biomet/81.4.673)\.
- \[13\]D\. Kohns and T\. Szendrei,Horseshoe prior Bayesian quantile regression, Journal of the Royal Statistical Society Series C: Applied Statistics, 73 \(2024\), pp\. 193–220,[https://doi\.org/10\.1093/jrsssc/qlad091](https://doi.org/10.1093/jrsssc/qlad091)\.
- \[14\]H\. Kozumi and G\. Kobayashi,Gibbs sampling methods for Bayesian quantile regression, Journal of Statistical Computation and Simulation, 81 \(2011\), pp\. 1565–1578,[https://doi\.org/10\.1080/00949655\.2010\.496117](https://doi.org/10.1080/00949655.2010.496117)\.
- \[15\]F\. Leibfried, V\. Dutordoir, S\. T\. John, and N\. Durrande,A tutorial on sparse Gaussian processes and variational inference, 2022,[https://doi\.org/10\.48550/arXiv\.2012\.13962](https://doi.org/10.48550/arXiv.2012.13962),[https://arxiv\.org/abs/2012\.13962](https://arxiv.org/abs/2012.13962)\.
- \[16\]N\. Meinshausen,Quantile regression forests, Journal of Machine Learning Research, 7 \(2006\), pp\. 983–999\.
- \[17\]H\. B\. Moss, S\. W\. Ober, and V\. Picheny,Inducing point allocation for sparse Gaussian processes in high\-throughput Bayesian optimisation, in Proceedings of The 26th International Conference on Artificial Intelligence and Statistics, PMLR, 2023, pp\. 5213–5230\.
- \[18\]K\. P\. Murphy,Probabilistic machine learning: an introduction, MIT Press, 2022\.
- \[19\]I\. Murray and R\. P\. Adams,Slice sampling covariance hyperparameters of latent Gaussian models, Advances in neural information processing systems, 23 \(2010\)\.
- \[20\]S\. W\. Ober, A\. Artemev, M\. Wagenländer, R\. Grobins, and M\. van der Wilk,Recommendations for baselines and benchmarking approximate Gaussian processes, 2024,[https://doi\.org/10\.48550/arXiv\.2402\.09849](https://doi.org/10.48550/arXiv.2402.09849),[https://arxiv\.org/abs/2402\.09849](https://arxiv.org/abs/2402.09849)\.
- \[21\]V\. Picheny, H\. Moss, L\. Torossian, and N\. Durrande,Bayesian quantile and expectile optimisation, in Proceedings of the Thirty\-Eighth Conference on Uncertainty in Artificial Intelligence, Eindhoven, Netherlands, 2022, PMLR, pp\. 1623–1633\.
- \[22\]M\. Plumlee and R\. Tuo,Building accurate emulators for stochastic simulations via quantile kriging, Technometrics, 56 \(2014\), pp\. 466–473,[https://doi\.org/10\.1080/00401706\.2013\.860919](https://doi.org/10.1080/00401706.2013.860919)\.
- \[23\]N\. Quadrianto, K\. Kersting, M\. D\. Reid, T\. S\. Caetano, and W\. L\. Buntine,Kernel conditional quantile estimation via reduction revisited, in 2009 Ninth IEEE International Conference on Data Mining, Miami Beach, Florida, USA, 2009, IEEE, pp\. 938–943,[https://doi\.org/10\.1109/ICDM\.2009\.82](https://doi.org/10.1109/ICDM.2009.82)\.
- \[24\]C\. E\. Rasmussen and C\. K\. I\. Williams,Gaussian processes for machine learning, Adaptive computation and machine learning, MIT Press, Cambridge, Massachusetts, USA, 3\. print ed\., 2006\.
- \[25\]C\. Sabater, O\. Le Maître, P\. M\. Congedo, and S\. Görtz,A Bayesian approach for quantile optimization problems with high\-dimensional uncertainty sources, Computer Methods in Applied Mechanics and Engineering, 376 \(2021\), p\. 113632,[https://doi\.org/10\.1016/j\.cma\.2020\.113632](https://doi.org/10.1016/j.cma.2020.113632)\.
- \[26\]J\. Sacks, W\. J\. Welch, T\. J\. Mitchell, and H\. P\. Wynn,Design and analysis of computer experiments, Statistical Science, 4 \(1989\), pp\. 409–423,[https://doi\.org/10\.1214/ss/1177012413](https://doi.org/10.1214/ss/1177012413)\.
- \[27\]E\. Snelson and Z\. Ghahramani,Sparse Gaussian processes using pseudo\-inputs, in Advances in Neural Information Processing Systems, vol\. 18, MIT Press, 2005\.
- \[28\]I\. Takeuchi, Q\. V\. Le, T\. D\. Sears, and A\. J\. Smola,Nonparametric quantile estimation, Journal of Machine Learning Research, 7 \(2006\)\.
- \[29\]P\. Thompson, Y\. Cai, R\. Moyeed, D\. Reeve, and J\. Stander,Bayesian nonparametric quantile regression using splines, Computational Statistics & Data Analysis, 54 \(2010\), pp\. 1138–1150,[https://doi\.org/10\.1016/j\.csda\.2009\.09\.004](https://doi.org/10.1016/j.csda.2009.09.004)\.
- \[30\]M\. Titsias,Variational learning of inducing variables in sparse Gaussian processes, in Proceedings of the Twelfth International Conference on Artificial Intelligence and Statistics, PMLR, 2009, pp\. 567–574\.
- \[31\]L\. Torossian, V\. Picheny, R\. Faivre, and A\. Garivier,A review on quantile regression for stochastic computer experiments, Reliability Engineering & System Safety, 201 \(2020\), p\. 106858,[https://doi\.org/10\.1016/j\.ress\.2020\.106858](https://doi.org/10.1016/j.ress.2020.106858)\.
- \[32\]S\. Vakili, H\. Moss, A\. Artemev, V\. Dutordoir, and V\. Picheny,Scalable Thompson sampling using sparse Gaussian process models, in Advances in Neural Information Processing Systems \(NeurIPS\), vol\. 34, Curran Associates, Inc\., 2021, pp\. 5631–5643\.
- \[33\]K\. Yu and R\. A\. Moyeed,Bayesian quantile regression, Statistics & Probability Letters, 54 \(2001\), pp\. 437–447,[https://doi\.org/10\.1016/S0167\-7152\(01\)00124\-9](https://doi.org/10.1016/S0167-7152(01)00124-9)\.
- \[34\]K\. Yu and J\. Zhang,A three\-parameter asymmetric Laplace distribution and its extension, Communications in Statistics \- Theory and Methods, 34 \(2005\), pp\. 1867–1879,[https://doi\.org/10\.1080/03610920500199018](https://doi.org/10.1080/03610920500199018)\.