A Review on Quantile Regression for Stochastic Computer Experiments

TL;DR

This paper empirically compares six quantile regression strategies for stochastic computer experiments, providing insights and guidelines for selecting the most suitable method based on problem characteristics.

Contribution

It introduces a comprehensive empirical evaluation of diverse quantile regression metamodels, offering practical guidelines for method selection in stochastic experiments.

Findings

Different metamodels perform variably depending on data characteristics.

Guidelines help choose the best quantile regression method for specific problem settings.

The study highlights the strengths and limitations of each approach.

Abstract

We report on an empirical study of the main strategies for quantile regression in the context of stochastic computer experiments. To ensure adequate diversity, six metamodels are presented, divided into three categories based on order statistics, functional approaches, and those of Bayesian inspiration. The metamodels are tested on several problems characterized by the size of the training set, the input dimension, the signal-to-noise ratio and the value of the probability density function at the targeted quantile. The metamodels studied reveal good contrasts in our set of experiments, enabling several patterns to be extracted. Based on our results, guidelines are proposed to allow users to select the best method for a given problem.

Figures

Tables

Table 1. Table 1: Summary of the metamodels, included the definition of the estimators, the associated numerical quantity and the related computation complexity. In this table L 𝐿 L and L η ∗ subscript 𝐿 superscript 𝜂 L_{\eta^{*}} are Hessian Lipschitz constants, M s subscript 𝑀 𝑠 M_{s} the number of multistarts (on the weights at hyperparameters fixed), d 𝑑 d the input dimension, N n ​ e ​ w subscript 𝑁 𝑛 𝑒 𝑤 N_{new} the size of the prediction set, n it subscript 𝑛 it n_{\operatorname*{it}} the number of iterations for the EM algorithm, k θ ( . , . ) k_{\theta}(.,.) the kernel function, B 𝐵 B a diagonal matrix, J 𝐽 J the total number of neurons, d ~ ~ 𝑑 \tilde{d} the number of covariables considered to find the best splitting point.

Method	Definition of ${\hat{q}}_{\tau }\left(x^{\ast }\right)$	Related quantity	Complexity
KN	$y_{\left(\left[K\tau \right]\right)}\left(x^{\ast }\right)$	The $K$-nearest points from $x^{\ast }$	$O\left(n{N}_{new}\left(d+\log n\right)\right)$
RF	$inf\{y_i:\hat{F}\left(y_i|X=x^{\ast }\right)\ge \tau \}$	$\hat{F}\left(y_i|X=x^{\ast }\right)=\sum _{i=1}^n{\overline{\omega }}_i\left(x^{\ast }\right){\mathrm{𝟙}}_{\left\{y_i\le y\right\}}$	$O\left(N\tilde{d}n{\log }^{2}n\right)$
	For a $3$ layer NN	With $w_i^{\left(h_1\right)}$ , $w_j^{\left(h_2\right)}$, $b_i^{\left(h_1\right)}$, $b_j^{\left(h_2\right)}$, $w^{\left(h_3\right)}$, $b^{\left(h_3\right)}$,
	$g_3(\sum _{j=1}^{J_2}{g}_2(\sum _{i=1}^{J_1}{g}_1(⟨w_i^{\left(h_1\right)},x^{\ast }⟩$	$1\le i\le J_1$, $1\le j\le J_2$, minimizing	$O\left(M_s{L}_{{\eta }^{\ast }}n{\left(Jd\right)}^3\right)$
NN	$+b_i^{\left(h_1\right)})w_j^{\left(h_2\right)}+b_j^{\left(h_2\right)})w^{\left(h_3\right)}+b^{\left(h_3\right)})$	$\frac{1}{n}\sum _{t=1}^n{l}_{\tau }\left(y_i-{\hat{q}}_{\tau }\left(x_i\right)\right)+\sum _{j,i}\frac{\lambda }{J_j}{\Vert w_i^{\left(h_j\right)}\Vert }^2$
		With $𝜶=({\alpha }_1,\dots ,{\alpha }_n)$ minimizing
		$\frac{1}{2}{𝜶}^T{K}_{x,x}^{\theta }𝜶-{𝜶}^T𝒚$
RK	$\sum _{i=1}^n{\alpha }_i{k}_{\theta }(x^{\ast },x_i)+b$	s.t $\frac{\tau -1}{\lambda n}\le {\alpha }_i\le \frac{\tau }{\lambda n},\forall 1\le i\le n$	$O\left(n^{7/2}\log \left(\frac{n}{ϵ}\right)\right)$
		and $\sum _{i=1}^n{\alpha }_i=0$ and $b$ the $\tau$-quantile of
		${\left(y_i-{\sum }_{j=1}^n{\alpha }_j{k}_{\theta }(x_i,x_j)\right)}_{1\le i\le n}$
		Maximizing the likelihood:
QK	$K_{x^{\ast },x}^{\theta }{\left(K_{x,x}^{\theta }+B\right)}^{-1}{\hat{q}}_{\tau }$	$p\left({𝒬}_{n^{\prime }}|{𝒳}_{n^{\prime }}\right)$	$O\left(L\left(d^3+n^{\prime 3}\right)+n^{\prime 3}\right)$
		${\hat{q}}_{\tau }\left(x_i\right)=q_{\tau }\left(x_i\right)+{\epsilon }_i{\epsilon }_i\sim 𝒩(0,B_{ii})$
		Approached solution that maximize:
VB	$K_{x^{\ast },x}^{\theta }{K}_{x,x}^{\theta -1}{\mu }_{\theta }$	$p\left({𝒴}_n|{𝒳}_n\right)$	$O\left(n_{it}\left(L\left(d^3+n^3\right)+n^3\right)\right)$
		$y_i=q_{\tau }\left(x_i\right)+\epsilon \epsilon \sim \text{ALP}(0,\sigma )$

Table 2. Table 2: Hyperparameters optimized on our benchmark.

Method	Hyperparameters
KN	number of neighbors $K\in {\mathbb{N}}^{\ast }$
RF	maximum size of the leaves $m_s\in {\mathbb{N}}^{\ast }$
NN	regularization $\lambda \in {ℝ}_{+}$, $J_1\in {\mathbb{N}}^{\ast }$
RK	regularization $\lambda \in {ℝ}_{+}$, lengthscales $\theta \in ℝ_{+}^{}{}_{}{}^d$
QK	length scale and variance $\theta \in ℝ_{+}^{}{}_{}{}^{d+1}$
VB	length scale and variance $\theta \in ℝ_{+}^{}{}_{}{}^{d+1}$

Table 3. Table 3: Summary of the characteristic of the problems.

		Test case 1	Test case 2	Test case 3	Test case 4	Sunflo
Dimension		$1$	2	1	9	9
Heteroscedasticity		very strong	very strong	very strong	strong	weak
Shape variation		very strong	weak	weak	weak	strong
pdf value	$\tau =0.1$	variable	small	globally very small	large	large
near the	$\tau =0.5$	variable	large	small	small	large
$\tau$-quantile	$\tau =0.9$	variable	very small	very large	very small	small

Table 4. Table 4: Data sizes for the different problems. The number in parentheses are the number of repetitions for QK.

Dimension	Data size (no repetitions)				Data size (with repetitions)
1	40	80	160	320	5 (8)	10 (8)	10 (16)	16 (20)
2	100	200	400	800	10 (10)	20 (10)	25 (16)	40 (20)
9	250	500	1000	2000	25 (10)	50 (10)	100 (10)	100 (20)

Equations

$$q_{\tau}=\min\big{\{}q\in\mathds{R}~{}:F(q)\geq\tau\big{\}}=\arg\min_{q\in\mathds{R}}\mathds{E}\big{[}l_{\tau}(Y-q)\big{]},$$

$$l_{\tau }(\xi )=(\tau -\text{mathds}{1}_{(\xi <0)})\xi ,\xi \in \text{mathds}R$$

$${\hat{q}}_{\tau }=\min \left\{y_i\in {\mathcal{Y}}_n:\hat{F}(y_i)\ge \tau \right\}$$

$${\hat{q}}_{\tau }=\arg \underset{q\in \text{mathds}R}{\min }\frac{1}{n}\sum _{i=1}^n{l}_{\tau }(y_i-q),$$

$$q_{\tau}(x)=\min\big{\{}q:F(q|X=x)\geq\tau\big{\}}=\arg\min_{q\in\mathds{R}}\mathds{E}\big{[}l_{\tau}(Y_{x}-q)\big{]},$$

$${\hat{q}}_{\tau }(x^{\ast })={\mathcal{Y}}_K^{x^{\ast }}([K\tau ]).$$

$$O(n{N}_{new}d)+O(n{N}_{new}\log n).$$

$$C(x_s)=\sum _{x_i\le x_s}(y_i-{\bar{Y}}_L{)}^2+\sum _{x_j>x_s}(y_j-{\bar{Y}}_R{)}^2,$$

$${\omega }_i(x^{\ast },t)=\frac{\text{mathds}{1}_{\{x_i\in \ell (x^{\ast },t)\}}}{\#\{j:x_j\in \ell (x^{\ast },t)\}},i=1,...,n$$

$${\bar{\omega }}_i(x^{\ast })=\frac{1}{N}\sum _{t=1}^N{\omega }_i(x^{\ast },t).$$

$$\hat{\mu }(x^{\ast })=\sum _{i=1}^n{\bar{\omega }}_i(x^{\ast })y_i.$$

$$\hat{F}(y|X=x^{\ast })=\sum _{i=1}^n{\bar{\omega }}_i(x^{\ast })\text{mathds}{1}_{\{y_i\le y\}}.$$

$$O\big{(}N\tilde{d}n\log^{2}(n)\big{)}$$

$$\hat{S}\in\arg\min_{s\in\mathcal{H}}\mathcal{R}_{e}[s]=\arg\min_{s\in\mathcal{H}}\dfrac{1}{n}\sum_{i=1}^{n}l\big{(}y_{i}-s(x_{i})\big{)}.$$

$${\mathcal{R}}_{r,e}[s]=\frac{1}{n}\sum _{i=1}^{n}l(y_i-s(x_i))+\lambda {\Vert s\Vert }^{\beta },$$

$$s(x^{*})=g_{3}\left(\sum_{j=1}^{J_{2}}g_{2}\Big{(}\sum_{i=1}^{J_{1}}g_{1}\big{(}\langle w_{i}^{(1)},x^{*}\rangle+b_{i}^{(1)}\big{)}w_{j,i}^{(2)}+b_{j}^{(2)}\Big{)}w_{1,j}^{(3)}+b^{(3)}\right).$$

$$\mathcal{R}_{r,e}[s]=\dfrac{1}{n}\sum_{i=1}^{n}l\big{(}y_{i}-s(x_{i})\big{)}+\lambda\sum_{j=1}^{H}\sum_{z=1}^{J_{j}}\left\|w_{z}^{(j)}\right\|^{2}.$$

$$l_{\tau }^{\eta }(\xi )=h^{\eta }(\xi )(\tau -\text{mathds}{1}_{\xi <0}),$$

$$h^{\eta }(\xi )=\{\begin{array}{r}\frac{{\xi }^2}{2\eta }\text{if}0\le |\xi |\le \eta \\ |\xi |-\frac{\eta }{2}\text{if}|\xi |\ge \eta .\end{array}$$

$${\mathbf{w}}_{\tau }^{\ast }=\underset{\mathbf{w}}{\mathrm{argmin}}\frac{1}{n}\sum _{t=1}^n{l}_{\tau }^{ϵ}(y_i-{\hat{q}}^{\mathbf{w}}(x_i))+\lambda \sum _{j=1}^H\sum _{i=1}^J{\Vert w_i^{(j)}\Vert }^2$$

$$\frac{{\mathcal{R}}_{r,e}^{\eta }({\mathbf{w}}_0)-{\mathcal{R}}_{r,e}^{\eta }({\mathbf{w}}^{\ast })}{\gamma }+{\log }_2{\log }_2({ϵ}_0/ϵ),$$

$$\frac{{\mathcal{R}}_{r,e}^{\eta }({\mathbf{w}}_0)-{\mathcal{R}}_{r,e}^{\eta }({\mathbf{w}}^{\ast })}{\gamma }+6.$$

$$O(M_s{L}_{{\eta }^{\ast }}n(Jd{)}^3),$$

$${\mathcal{R}}_{r,e}[s]=\frac{1}{n}\sum _{i=1}^{n}l(y_i-s(x_i))+\frac{\lambda }{2}{\Vert s\Vert }_{\mathcal{H}}^2,$$

$$k(x,x^{\prime })={⟨\Phi (x^{\prime }),\Phi (x)⟩}_{\mathcal{H}}.$$

$$s(x)={⟨s,k(.,x)⟩}_{\mathcal{H}}$$

$${\mathcal{H}}_{|X}^{\theta }=\mathrm{span}\{\Phi (x_i):i=1,...,n\}\text{with}{\Vert s\Vert }_{{\mathcal{H}}_{|X}}^2=\sum _{i=1}^n\sum _{j=1}^n{\alpha }_i{\alpha }_j{k}_{\theta }(x_j,x_i).$$

$$\hat{S}(x)=\sum _{i=1}^n{\alpha }_i{k}_{\theta }(x,x_i).$$

$$\dfrac{1}{n}\sum_{i=1}^{n}l\Big{(}y_{i}-\big{(}\sum_{j=1}^{n}\alpha_{j}k_{\theta}(x_{i},x_{j})\big{)}\Big{)}+\dfrac{\lambda}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\alpha_{i}\alpha_{j}k_{\theta}(x_{j},x_{i}).$$

$$q_{\tau }(x)=g(x)+b\text{with}g\in {\mathcal{H}}_{\theta }\text{with}b\in \text{mathds}R,$$

Full text

1 Introduction

Computer simulation models are now essential for performance evaluation, quality control and uncertainty quantification to assess decisions in complex systems. These computer simulators generally model systems depending on multiple input variables that can be divided into two categories: the controllable variables and the uncontrollable variables.

For example, in pharmacology, the optimal drug dosage depends on the drug formulation but also on the targeted individual (genetics, age, sex) and environmental interactions. The shelf life and performance of a manufacturing device depend on its design, but also on its environment and on some uncertainties due to the manufacturing process. The plant growth and yield depend on the genes of the plant and on the gardening techniques but also on the weather and potential diseases.

Evaluating the influence of the controllable and uncontrollable variables directly on the real-life problems can be costly and tedious. One solution is to encapsulate the systems into a computer simulation model which would reduce the cost and the time required for each test (see [40] and reference therein for computer experiments applied to clinical trials, see [35] for a computer simulation model applied to industrial design and [91, 19] for computer experiments applied to crop production).

In such computer simulation models, the links between the inputs and outputs may be too complex to be fully understood or to be formulated in a closed form. In this case, the system can be considered as a black box and formalized by an unknown function: $\Psi:\mathcal{X}\times\Omega\rightarrow\mathds{R}$, where $\mathcal{X}\subset\mathds{R}^{d}$ denotes the compact space of controllable variables, and $\Omega$ denotes a probability space representing the uncontrollable variables. Note that some black boxes have their outputs in a multidimensional space but this aspect is beyond the scope of the present paper.

Based on the standard constraints encountered in computer experiments, throughout this paper, we assume that the function is only accessible through pointwise evaluations $\Psi(x,\omega)$; no structural information is available regarding $\Psi$; in addition, we take into account the fact that evaluations may be expensive, which drastically limits the number of possible calls to $\Psi$.

If the space $\Omega$ is small enough (or highly structured), it may be possible to work directly on the space $\mathcal{X}\times\Omega$ (see [42] for example). Often, however, $\Omega$ is too complex and working on the joint space is intractable. This is the case for intrinsic stochastic simulators (see [48, 54] for examples of biological systems, where the stochasticity is driven by stochastic equations such as the Fokker-Planck or the chemical Langevin equations), or for simulators associated with a very large space $\Omega$ (see [19] for the crop model SUNFLO that considers $5$ weather indicators on $130$ days, $i.e$ $\Omega$ is a space of dimension $650$).

In this paper we consider the case where $\Omega$ is too complex. We assume that $\Omega$ has any structure and its contribution is considered as random. In contrast to deterministic systems, for any fixed $x$, $\Psi(x,\cdot)$ is considered as a random variable of distribution $\mathds{P}_{x}$; hence, such systems are often referred to as stochastic black boxes. In order to understand the behavior of the system of interest or to take optimal decisions, information is needed about the output distribution. An intuitive approach is to use a simple Monte-Carlo technique and to evaluate $\Psi(x,\omega_{1}),\ldots,\Psi(x,\omega_{n})$ to extract statistical moments, the empirical cumulative distribution function, etc. Unfortunately, such a stratified approach is not efficient when evaluating $\Psi$ is expensive.

Instead, we focus on surrogate models (also referred to as metamodels or statistical emulators), which are appropriate approaches in small data settings [84, 94]. Among the vast choice of surrogate models, the most popular ones include regression trees, Gaussian processes, support vector machines and neural networks. In the framework of stochastic black boxes, the standard approach consists in estimating the conditional expectation of $\Psi$. This case has been extensively treated in the literature and many applications, including Bayesian optimization [77], have been developed. However the conditional expectation is risk-neutral, as it does not provide information about the variability of the output, whereas pharmacologists, manufacturers, asset managers, data scientists and agronomists may need to quantify potential losses associated with their decisions.

Risk information can be introduced via heteroscedastic Gaussian processes [43, 47] in which the unknown distribution is estimated by a surrogate expectation-variance model. However, such approaches usually imply that the shape of the distribution (e.g. normal) is the same for all $x\in\mathcal{X}$. Another possible approach would be to learn the unknown distribution with no strong structural hypotheses [61, 37, 31], but this requires a large number of evaluations of $\Psi$. Here, we focus on the conditional quantile estimation of order $\tau$, a flexible way to tackle cases in which the distribution of $\Psi(x,.)$ varies markedly in spread and shape with respect to $x\in\mathcal{X}$, and a classical risk-aware tool in decision theory [71].

1.1 Paper Overview

Many metamodels originally designed to estimate conditional expectations have been adapted to estimate the conditional quantile. However, despite extensive literature on estimating the quantile in the presence of spatial structure, few studies have reported on the constraints associated with stochastic black boxes. The performance of a metamodel with high dimension input is treated in insufficient details, performance based on the number of points has rarely been tackled and, to our knowledge, dependence on specific aspects of the quantile functions has never been studied. The aim of the present paper is to review quantile regression methods under standard constraints related to the stochastic black box framework, so as to provide information on the performance of the selected methods, and to recommend which metamodel to use depending on the characteristics of the computer simulation model and the data.

A comprehensive review of quantile regression is of course beyond the scope of the present work. We limit our review to the approaches that are best suited for our framework, while ensuring the necessary diversity of metamodels. In particular, we have chosen six metamodels that are representative of three main categories: approaches based on statistical order (K-nearest neighbors [KN] regression [7] and random forest [RF] regression [58]), functional or frequentist approaches (neural networks [NN] regression [18] and regression in reproducing kernel Hilbert space [RK] [86]), and Bayesian approaches based on Gaussian processes (Quantile Kriging [QK] [64] and the variational Bayesian [VB] regression [1]). Each category has some specificities in terms of theoretical basis, implementation and complexity. We begin this presentation by describing the methods in full in sections 3, 4 and 5.

In order to identify the relevant areas of expertise of the different metamodels, an original benchmark system is designed based on four toy functions and an agronomical model [19]. The dimension of the problems ranges from $1$ to $9$ and the number of observations from $40$ to $2000$. Particular attention is paid to the performance of each metamodel according to the size of the learning set, the value of the probability density function at the targeted quantile $f(.,q_{\tau})$ and the dimension of the problem. Sections 6 and 7 describe the benchmark system and detail its implementation, with particular focus on the tuning of the hyperparameters of each method. Full results and discussion are to be found in Sections 8 and 9, respectively.

2 Quantile emulators and design of experiments

We first provide the necessary definitions, objects and properties related to the quantile. The quantile of order $\tau\in(0,1)$ of a random variable $Y$ can be defined either as the (generalized) inverse of a cumulative distribution function (CDF), or as the solution to an optimization problem:

[TABLE]

$F(\cdot)$ is the CDF of $Y$ and

[TABLE]

is the so-called pinball loss [45] (Figure 1). In the following, we only consider situations in which $F$ is strictly increasing and continuous.

Given a finite observation set $\mathcal{Y}_{n}=(y_{1},...,y_{n})$ composed of i.i.d samples of $Y$, the empirical estimator of $q_{\tau}$ can thus be introduced in two different ways:

[TABLE]

or

[TABLE]

where $\hat{F}$ denotes an estimator of the CDF function. In (3) $\hat{q}_{\tau}$ coincides with an order statistic. For example, if $\hat{F}$ is the empirical CDF function then $\hat{q}_{\tau}=y_{([n\tau])}$, where $[n\tau]$ represents the smallest integer greater than or equal to $n\tau$ and $y_{(k)}=\mathcal{Y}_{n}(k)$ is the $k$-th smallest value in the sample $\{y_{1},\ldots,y_{n}\}$. The estimators (4) and (3) may coincide, but are in general not equivalent.

Similarly to (1), the conditional quantile of order $\tau\in(0,1)$ can be defined in two equivalent ways:

[TABLE]

where $Y_{x}$ is a random variable of distribution $\mathds{P}_{x}$ and $F(.|X=x)$ is the CDF of $Y_{x}$.

In a quantile regression context, one only has access to a finite observation set \mathcal{D}_{n}=\big{\{}(x_{1},y_{1}),\allowbreak\ldots,\allowbreak(x_{n},y_{n})\big{\}}=(\mathcal{X}_{n},\mathcal{Y}_{n}) with $\mathcal{X}_{n}$ a $n\times d$ matrix. Estimators for (5) are either based on the order statistic as in (3) (section 3), or on a minimizer of the pinball loss as in (4) (sections 4 and 5). Throughout this work, the observation set $\mathcal{D}_{n}$ is fixed (we do not consider a dynamic or sequential framework). Following the standard approach used in computer experiments, the training points $x_{i}$ are chosen according to a space-filling design [20] over a hyperrectangle. In particular, we assume that there are no repeated experiments: $x_{i}\neq x_{j}$, $\forall i\neq j$; most of the methods chosen in this survey (KN, RF, RK, NN, VB) work under that setting.

However, as a baseline approach, one may decide to use a stratified experimental design $\mathcal{D}_{n^{\prime},r}$ with $r$ i.i.d samples for a given $x_{i}$, $i=1,..,n^{\prime}$, extract pointwise quantile estimates using (3) and fit a standard metamodel to these estimates. The immediate drawback is that for the same budget ($n^{\prime}\times r=n$) such experimental designs cover much less of the design space than a design with no repetition. The QK method is based on this approach.

3 Methods based on order statistics

A simple way to compute a quantile estimate is to take an order statistic of an i.i.d. sample. A possible approach is to emulate such a sample by selecting all the data points in the neighborhood of the query point $x$, and then by taking the order statistic of this subsample as an estimator for the conditional quantile. One may simply choose a subsample of $\mathcal{D}_{n}$ based on a distance defined on $\mathcal{X}$: this is what the $K$-nearest neighbors approach does. It is common to use KN based on the Euclidean distance but of course any other distance can be used, such as Mahalanobis [93] or weighted Euclidean distance [30]. Alternatively, one may define a notion of neighborhood using some space partitioning of $\mathcal{X}$. That includes all the decision tree methods [16], in particular regression trees, bagging or random forest [58].

3.1 $K$-nearest neighbors

The $K$-nearest neighbors method was first proposed for the estimation of conditional expectations [82, 83]. Its extension to the conditional quantile estimation can be found in [7].

3.1.1 Quantile regression implementation

Define $\mathcal{X}_{\operatorname*{test}}$ as the set of query points. KN works as follows: for each $x^{*}\in\mathcal{X}_{\operatorname*{test}}$ define $\mathcal{X}^{K}(x^{*})$ the subset of $\mathcal{X}_{n}$ containing the $K$ points that are the closest to the query point $x^{*}$. Define $\mathcal{Y}_{K}^{x^{*}}$ the associated outputs, and define $\hat{F}^{K}(y|X=x^{*})$ as the associated empirical CDF. Following (3), the conditional quantile of order $\tau$ can be defined as the statistical order

[TABLE]

Algorithm 1 details the implementation of the KN method.

3.1.2 Computational complexity

For a naive implementation of such an estimator, one needs to compute $n\times N_{new}$ distances, where $N_{new}$ is the number of query points, hence for a cost in $O(nN_{new}d)$. Moreover, sorting $n$ distances in order to extract the $K$ nearest points has a cost in $O(nN_{new}\log n)$. Combining the two operations implies a complexity of order

[TABLE]

Note that some algorithms have been proposed in order to reduce the computational time, for example by using GPUs [34] or by using tree search algorithms [4].

3.2 Random forests

Random forests were introduced by Breiman [15] for the estimation of conditional expectations. They have been used successfully for classification and regression, especially with problems where the number of variables is much larger than the number of observations [27].

3.2.1 Overview

The basic element of random forests is the regression tree $T$, a simple regressor built via a binary recursive partitioning process. Starting with all data in the same partition $i.e$ $\mathcal{X}$, the following sequential process is applied. At each step, the data is split into two, so that $\mathcal{X}$ is partitioned in a way that it can be represented by a tree as it is presented Figure 2.

Several splitting criteria can be chosen (see [41]). In [58], the splitting point $x_{S}$ is the data point that minimizes

[TABLE]

where $\bar{Y}_{L}$ and $\bar{Y}_{R}$ are the mean of the left and right sub-populations, respectively. Equation (7) applies when the $x$’s are real-valued. In the multidimensional case, the dimension $d_{S}$ in which the split is performed has to be selected. The split then goes through $x_{S}$ and perpendicularly to the direction $d_{S}$. There are several rules to stop the expansion of $T$. For instance, the process can be stopped when the population of each cell is inferior to a minimal size $nodesize$: then, each node becomes a terminal node or leaf. The result of the process is a partition of the input space into hyperrectangles $R(T)$. Like the KN method, the tree-based estimator is constant on each neighborhood. The hope is that the regression trees automatically build neighborhoods from the data that should be adapted to each problem.

Despite their simplicity of construction and interpretation, regression trees are known to suffer from a certain rigidity and a high variance (see [14] for more details). To overcome this drawback, regression trees can be used with ensemble methods like bagging. Instead of using only one tree, bagging creates a set of tree $\mathcal{T}_{N}=\left\{T^{1},..,T^{N}\right\}$ based on a bootstrap version \mathcal{D}_{N,n}=\big{\{}\big{(}(x_{1_{t}},y_{1_{t}}),...,(x_{n_{t}},y_{n_{t}})\big{)}\big{\}}_{t=1}^{N} of $\mathcal{D}_{n}$. Then the final model is created by averaging the results among all the trees.

Bagging reduces the variance of the predictor, as the splitting criterion has to be optimized over all the input dimensions, but computing (7) for each possible split is costly when the dimension is large. The random forest algorithm, a variant of bagging, constructs an ensemble of weak learners based on $\mathcal{D}_{N,n}$ and aggregates them. Unlike plain bagging, at each node evaluation, the algorithm uses only a subset of $\tilde{d}$ covariables for the choice of the split dimensions. Because the $\tilde{d}$ covariables are randomly chosen, the result of the process is a random partition $R(t)$ of $\mathcal{X}$ constructed by the random tree $T^{t}$.

3.2.2 Quantile prediction

We present the extension proposed in [58] for conditional quantile regression. Let us define $\ell(x^{*},t)$ the leaf obtained from the tree $t$ containing a query point $x^{*}$ and

[TABLE]

[TABLE]

The $\bar{\omega}_{i}(x^{*})$’s represent the weights illustrating the “proximity” between $x^{*}$ and $x_{i}$. In the classical regression case, the estimator of the expectation is:

[TABLE]

In [58] the conditional quantile of order $\tau$ is defined as in (3) with the CDF estimator defined as

[TABLE]

Algorithm 2 details the implementation of the RF method.

3.2.3 Computational complexity

Assuming that the value of (7) can be computed sequentially for consecutive thresholds, the RF computation burden lies in the search of the splitting point that implies sorting the data. Sorting $n$ variables has a complexity in $O(n\log n)$. Thus, at each node the algorithm finds the best splitting points considering only $\tilde{d}\leq d$ covariables (classically $\tilde{d}=d/3$). This implies a complexity of O\big{(}\tilde{d}n\log n\big{)} per node. In addition, the depth of a tree is generally upper bounded by $\log n$. Then the computational cost of building a forest containing $N$ trees under the criterion (7) is

[TABLE]

[53, 95]. One may observe that RF are easy to parallelize and that contrary to KN the prediction time is very small once the forest is built.

4 Approaches based on functional analysis

Functional methods search directly for the function mapping the input to the output in a space fixed beforehand by the user. With this framework, estimating any functional $S$ of the conditional distribution implies selecting a loss $l$ (associated to $S$) and a function space $\mathcal{H}$. Thus, the estimator $\hat{S}\in\mathcal{H}$ is obtained as the minimizer of the empirical risk $\mathcal{R}_{e}$ associated to $l$, $i.e.$

[TABLE]

The functional space $\mathcal{H}$ must be chosen flexible enough to extract some signal from the data. In addition, $\mathcal{H}$ needs to have enough structure to make the optimization procedure feasible (at least numerically). In the literature, several formalisms such as linear regression [76], spline regression [56], support vector machine [92], neural networks [9] or deep neural networks [74] use structured functional spaces with different levels of flexibility.

However, using a too large $\mathcal{H}$ can lead to overfitting, $i.e.$ return predictors that are good only on the training set and generalize poorly. Overcoming overfitting requires some regularization (see [75, 100, 101] for instance), defining for example the regularized risk

[TABLE]

where $\lambda\in\mathds{R}^{+}$ is a penalization factor, $\beta\in\mathds{R}^{+}$ and $\left\|.\right\|$ is either a norm for some methods (Section 4.2) or a measure of variability for others (Section 4.1). The parameter $\lambda$ plays a major role, as it allows to tune the balance between bias and variance.

Classically, squared loss is used: it is perfectly suited to the estimation of the conditional expectation. Using the pinball loss (Eq. 2) instead allows to estimate quantiles. In this section we present two approaches based on Equation (11) with the pinball loss. The first one is regression using artificial neural networks (NN), a rich and versatile class of functions that has shown a high efficiency in several fields. The second approach is the generalized linear regression in reproducing kernel Hilbert spaces (RK). RK is a non-parametric regression method that has been much studied in the last decades (see [80]) since it appeared in the core of learning theory in the 1990’s [75, 92].

4.1 Neural Networks

Artificial neural networks have been successfully used for a large variety of tasks such as classification, computer vision, music generation, and regression [9]. In the regression setting, feed-forward neural networks have shown outstanding achievements. Here we present quantile regression neural network [18] which is an adaptation of the traditional feed-forward neural network.

4.1.1 Overview

A feed-forward neural network is defined by its number of hidden layers $H$, its numbers of neurons per layer $J_{h},1\leq h\leq H$, and its activation functions $g_{h}$, $h=1,\dots,H$. Given an input vector $x\in\mathds{R}^{d}$ the information is fed to the hidden layer $1$ composed of a fixed number of neurons $J_{1}$. For each neuron $N_{i}^{(1)}$, $i=1,..,J_{1}$, a scalar product (noted $\left\langle.~{},.\right\rangle$) is computed between the input vector $x=(x_{1},...,x_{d})\in\mathds{R}^{d}$ and the weights $w_{i}^{(1)}=(w_{i,1}^{(1)},...,w_{i,d}^{(1)})\in\mathds{R}^{d}$ of the $N_{i}^{(1)}$ neurons. Then a bias term $b_{i}^{(1)}\in\mathds{R}$ is added to the result of the scalar product. The result is composed with the activation function $g_{1}$ (linear or non-linear) which is typically the sigmoid or the ReLu function [74] and the result is given to the next layer where the same operation is processed until the information comes out from the outpout layer. For example, the output of a $3$-layers NN at $x^{*}$ is given by

[TABLE]

The corresponding architecture can be found in Figure 3.

The architecture of the NN defines $\mathcal{H}$. Finding the right architecture is a very difficult problem which will not be treated in this paper. However a classical implementation procedure consists of creating a network large enough (able to overfit) and then using techniques such as early stopping, dropout, bootstrapping or risk regularization to avoid overfitting [79]. In [18], the following regularized risk is used:

[TABLE]

4.1.2 Quantile regression

Minimizing Equation (13) (with respect to all the weights and biases) is in general challenging, as $\mathcal{R}_{r,e}$ is a highly multimodal function. It is mostly tackled using derivative-based algorithms and multi-starts ($i.e$ launching the optimization procedure $M_{s}$ times with different starting points). In the case of quantile estimation, the loss function is non-differentiable at the origin, which may cause problems to some numerical optimization procedures. To address this issue, [18] introduced a smooth version of the pinball loss function, defined as:

[TABLE]

where

[TABLE]

Note that if the optimizer is based on a first order method such as [44], then the transfer function does not require continuous derivatives. But using a second order method as it is done in the original paper implies the loss function to be twice differentiable with respect to the weights of the neural network. Then transfer functions such as logistic or hyperbolic tangent functions should be used over piecewise linear ones such as the ReLu or the PReLU functions [66].

Let us define $\boldsymbol{w}$ the list containing the weights and bias of the network. To find $\boldsymbol{w}^{*}$, a minimizer of $\mathcal{R}_{r,e}$, the idea is to solve a series of problems using the smoothed loss instead of the pinball one with a sequence $E_{K}$ corresponding to $K$ decreasing values of $\eta$. The process begins with the optimization with the larger value $\eta_{1}$. Once the optimization converges, the optimal weights are used as the initialization for the optimization with $\eta_{2}$, and so on. The process stops when the weights based on $l_{\tau}^{\eta_{K}}$ are obtained. Finally, $\hat{q}_{\tau}(x^{*})$ is given by the evaluation of the optimal network at $x^{*}$. Algorithm 3 details the implementation of the NN method.

4.1.3 Computational complexity

In [18] the optimization is based on a Newton method. Thus, the procedure needs to inverse a Hessian matrix. Without sophistication, its cost is $O(s_{\operatorname*{pb}}^{3})$ with $s_{\operatorname*{pb}}$ the size of the problem $i.e$ the number of parameters to optimize. Note that using a high order method makes sense here because NN has few parameters (in contrast to deep learning methods). Moreover providing an upper bound on the number of iterations needed to reach an optimal point may be really hard in practice because of the non convexity of (13). In the non-convex case, there is no optimality guaranty and the optimization could be stuck in a local minima. However, it can be shown that the convergence near a local optimal point is at least super linear (see [12] Eq. (9.33)) and may be quadratic (if the gradient is small). It implies, for each $\eta$, the number of iterations until $\mathcal{R}_{r,e}^{\eta}(\boldsymbol{w})-\mathcal{R}_{r,e}^{\eta}(\boldsymbol{w}^{*})\leq\epsilon$ is bounded above by

[TABLE]

with $\gamma$ the minimal decreasing rate, $\epsilon_{0}=2M_{\eta}^{3}/L_{\eta}^{2}$, $M_{\eta}$ the strong convexity constant of $\mathcal{R}_{r,e}^{\eta}$ near $\boldsymbol{w}^{*}$ and $L_{\eta}$ the Hessian Lipschitz constant (see [12] page 489). As $\log_{2}\log_{2}(\epsilon_{0}/\epsilon)$ increases very slowly with respect to $\epsilon$, it is possible to bound the number of iterations $N$ typically by

[TABLE]

That means, near an optimal point, the complexity is $O(L_{\eta}n(Jd)^{3})$, with $J$ the total number of neurons. Then using a multistart procedure implies a complexity of

[TABLE]

with $L_{\eta^{*}}=\max_{\eta_{1},\dots,\eta_{K}}L_{\eta}$ .

4.2 Regression in RKHS

Regression in RKHS was introduced for classification via Support Vector Machine by [25, 39], and has been naturally extended for the estimation of the conditional expectation [29, 70]. Since, many applications have been developed (see [80, 75] for some examples), here we present the quantile regression in RKHS [86, 73].

4.2.1 RKHS introduction and formalism

Under the linear regression framework, $S$ is assumed to be under the form $S(x)=x^{T}\alpha$, with $\alpha$ in $\mathds{R}^{d}$. To stay in the same vein while creating non-linear responses, one can map the input space $\mathcal{X}$ to a space of higher dimension $\mathcal{H}$ (named the feature space), thanks to a feature map $\Phi$. Nevertheless, although large dimensional spaces enable more flexibility, it may introduce difficulties during the estimation procedure. A good trade-off can be found by combining regularized empirical risk

[TABLE]

and the RKHS formalism that is based on the representer theorem and the so-called kernel trick. Define $\mathcal{H}$ as a $\mathds{R}$-Hilbert functional space and let us assume there exists a symmetric definite positive function $k:\mathcal{X}\times\mathcal{X}\rightarrow\mathds{R}$ such that:

[TABLE]

Under this setting, $\mathcal{H}$ is a RKHS with the reproducing kernel $k$, that means $\Phi(x)=k(.,x)\in\mathcal{H}$ for all $x\in\mathcal{X}$ and the reproducing property

[TABLE]

holds for all $s\in\mathcal{H}$ and all $x\in\mathcal{X}$. It can be shown that working with a fixed kernel $k$ is equivalent to work with its associated functional Hilbert space. Note that the kernel choice is based on kernel properties or assumptions made on the functional space. See for instance [80], chapter 4, for some kernel definitions and properties. In the following, $\mathcal{H}_{\theta}$ and $k_{\theta}$ denote respectively a RKHS and its kernel associated to the hyperparameters vector $\theta$. $K_{x,x}^{\theta}\in\mathcal{R}^{n\times n}$ is the kernel matrix obtained via $K_{x,x}^{\theta}(i,j)=k_{\theta}(x_{i},x_{j})$.

From a theoretical point of view, the representer theorem implies that the minimizer $\hat{S}$ of (15) lives in

[TABLE]

Combining this result to the definition (16) and the reproducing property (17), it is possible to write $\hat{S}$ as:

[TABLE]

Hence, the estimation procedure becomes an optimization problem over $n$ coefficients $\boldsymbol{\alpha}=(\alpha_{1},\alpha_{2},...,\alpha_{n})\in\mathds{R}^{n}$. More precisely, finding $\hat{S}$ is equivalent to minimize with respect to $\boldsymbol{\alpha}$ the quantity

[TABLE]

4.2.2 Quantile regression

Quantile regression in RKHS was introduced by [86], followed by several authors [50, 81, 23, 24, 73]. Quantile regression has two specificities compared to the general case. Firstly the loss $l$ is defined as the pinball. Secondly, to ensure the quantile property, the intercept is not regularized. More precisely, we assume that

[TABLE]

and we consider the empirical regularized risk

[TABLE]

Thus, the representer theorem implies that $\hat{q}_{\tau}$ can be written under the form

[TABLE]

for a new query point $x^{*}$. Since (19) cannot be minimized analytically, a numerical minimization procedure is used. [25] followed by [86] introduced non-negative variables $\xi^{(*)}\in\mathds{R}^{+}$ to transform the original problem into

[TABLE]

subject to

[TABLE]

and

[TABLE]

Using a Lagrangian formulation, it can be shown (see [80] for instance) that minimizing $\mathcal{R}_{r,e}$ is equivalent to the problem:

[TABLE]

It is a quadratic optimization problem under linear constraint, for which many efficient solvers exist.

The value of $b$ may be obtained from the Karush-Kuhn-Tucker slackness condition or fixed independently of the problem. A simple way to do so is to choose $b$ as the $\tau$-quantile of \big{(}y_{i}-\sum_{j=1}^{n}\alpha_{j}k_{\theta}(x_{i},x_{j})\big{)}_{1\leq i\leq n}. Algorithm 4 details the implementation of the RK method.

4.2.3 Computational complexity

Let us notice two things. Firstly, the minimal upper bound complexity for solving (20) is $O(n^{3})$. Indeed solving (20) without the constraints is easier and it needs $O(n^{3})$. Secondly the optimization problem (20) is convex, thus the optimum is global.

There are two main approaches for solving (20), the interior point method [12] and the iterative methods like libSVM [21]. The interior point method is based on the Newton algorithm, one method is the barrier method (see [12] page 590). It is shown that the number of iterations $N$ for reaching a solution with precision $\epsilon$ is $O(\sqrt{n}\log(\frac{n}{\epsilon})).$ Moreover each iteration of a Newton type algorithm costs $O(n^{3})$ because it needs to inverse a Hessian. Thus, the complexity of an interior point method for finding a global solution with precision $\epsilon$ is

[TABLE]

On another hand, iterative methods like libSVM transform the main problem into a smaller one. At each iteration the algorithm solves a subproblem in $O(n)$. Contrary to the interior point methods, the number of iterations depends explicitly on the matrix $K_{x,x}^{\theta}$. [51] shows that the number of iterations is

[TABLE]

where $\kappa(K_{x,x}^{\theta})=\lambda_{\max}(K_{x,x}^{\theta})/\lambda_{\min}(K_{x,x}^{\theta})$. Note that $\kappa(K_{x,x}^{\theta})$ depends on the type of the kernel, it evolves in $O(n^{s^{\prime}})$ with $s^{\prime}>1$ an increasing value of the regularity of $k_{\theta}$ [22]. For more information about the eigenvalues of $K_{x,x}^{\theta}$ one can consult [13].

To summarize, it implies that the complexity of the libSVM method has an upper bound higher than the interior point algorithm. However, these algorithms are known to converge pretty fast. In practice, the upper bound is almost never reached, and thus the most important factor is the cost per iteration, rather than the number of iterations needed. This is the reason why libSVM is popular in this setting.

5 Bayesian approaches

Bayesian formalism has been used for a wide class of problems such as classification and regression [67], model averaging [11] and model selection [65]. The first Bayesian quantile regression model was introduced in [97] where the authors worked under linear frameworks combined with improper uniform priors. Nevertheless the linear hypothesis may be too restrictive to treat the stochastic black box setting. As an alternative [85] introduced a mixture modeling framework called Dirichlet process to perform nonlinear quantile regression. However the inference is performed with MCMC methods (see [36, 33] for instance), a procedure that is often costly. A possibility to overcome that drawback is the use of Gaussian process (GP). GPs are powerful in a Bayesian context because of their flexibility and their tractability (GPs are only characterized by their mean $m$ and covariance $k_{\theta}$). Using GPs, a possible approach is to use a joint modeling of the mean and variance, assuming that the distribution is Gaussian everywhere, and then to extract the quantiles of interest as it is done in [43, 47]. Nevertheless this strategy may introduce a high bias when the true output distribution is not Gaussian.

In this section we present QK and VB, two approaches that use GP as a prior for $q_{\tau}$. Contrary to the joint modelling approach, here the GP prior for $q_{\tau}$ does not imply any structure on the output distribution, which allows the creation of very flexible quantile models.

5.1 Kriging

Kriging is based on the hypothesis that

[TABLE]

Which means for every finite set $(x_{1},\cdots,x_{T})$, the output \big{(}S(x_{1}),\cdots,S(x_{T})\big{)} is multivariate Gaussian. Here $m$ is the mean of the process and $k_{\theta}$ is a kernel function also known as the covariance function. Note that in the sequel we take $m=0$ in order to simplify the computations and notations. The covariance function conveys many properties of the process, so its choice should depend on the assumptions made on $S$. A first classical assumption is that the GP is stationary, $i.e$ the correlation between two inputs does not depends on the location but only on the distance between the points. Then, different class of stationary kernels produce GPs with different regularities. The class of Matérn kernels is very convenient because it depends on a regularity hyperparameter that enables the end user to adapt its prior to his regularity assumptions (see [67] for more details). For example, the Matérn $5/2$ kernel is defined as

[TABLE]

where $\rho>0$ and

[TABLE]

with $\Lambda_{\theta}$ a diagonal matrix with diagonal terms the inverses of the $d$ squared length scales $\theta_{i}$, $i=1,\cdots,d$.

In addition to the assumption 21, let us assume $y_{i}$ is observed with noise such that

[TABLE]

with $\varepsilon_{i}\sim\mathcal{N}(0,\sigma_{i}^{2})$. As a consequence the associated likelihood is Gaussian, $i.e.$

[TABLE]

with $\boldsymbol{\sigma}=(\sigma_{1},\cdots,\sigma_{n}).$ Because \big{(}\mathcal{Y}_{n},S(x_{*})\big{)}^{T} is a Gaussian vector of zero mean and covariance

[TABLE]

the distribution of $S(x_{*})$ knowing $\mathcal{Y}_{n}$ is still Gaussian (see [67] appendix A.2 for more details). Thus, it is possible to provide the distribution a posteriori for Kriging regression as

[TABLE]

with

[TABLE]

[TABLE]

As in Section 4.2, the covariance functions are usually chosen among a set of predefined ones (for example the Matérn $5/2$ see Eq.22), that depend on a set of hyperparameters $\boldsymbol{\theta}\in\mathds{R}^{d+1}$. The best hyperparameter $\boldsymbol{\theta}^{*}$ can be selected as the maximizer of the marginal likelihood. More precisely it follows (see [67] for instance) that

[TABLE]

where $|K|$ is the determinant of the matrix $K$. Maximizing this likelihood with respect to $\boldsymbol{\theta}$ is usually done using derivative-based algorithms, although the problem is non-convex and known to have several local maxima.

Different estimators of $S$ may be extracted based on (26). Here $\hat{S}$ is fixed as $\bar{S}_{\theta^{*}}$. Note that this classical choice is made because the maximum a posteriori of a Gaussian distribution coincides with its mean.

5.1.1 Quantile kriging

As $q_{\tau}$ is a latent quantity, the solution proposed in [64] is to consider the sample $\mathcal{D}_{n^{\prime},r}$ corresponding to a design of experiments with $n^{\prime}$ different points that are repeated $r$ times in order to obtain quantile observations. For each $x_{i}\in\mathcal{X}$, $1\leq i\leq n^{\prime}$, let us define:

[TABLE]

and

[TABLE]

Following [64], let us assume that

[TABLE]

Note that from a statistical point of view Assumption (28) is wrong because the distribution is asymetric around the quantile but asymptotically consistent as illustrated by the central limit theorem for sample quantiles. The resulting estimator is

[TABLE]

with \mathcal{Q}_{n^{\prime}}=\big{(}\hat{q}_{\tau}(x_{1}),\dotsb,\hat{q}_{\tau}(x_{n^{\prime}})\big{)} and $B=\text{diag}(\sigma_{1}^{2},\ldots,\sigma_{n}^{\prime 2})$.

There are several possibilities to evaluate the noise variances $\sigma_{i}^{2}$. Here we choose to use a bootstrap technique (that is, generate bootstrapped samples of $y_{i,r}$, compute the corresponding $\hat{q}_{\tau}(x_{i})$ values and take the variance over those values as the noise variance) but it is possible to use the central limit theorem as it is presented in [5]. The hyperparameters are selected based on $(\ref{likel})$ changing $\mathcal{Y}_{n}$ by $\mathcal{Q}_{n^{\prime}}$. Algorithm 5 details the implementation of the QK method.

5.1.2 Computational complexity

Optimizing (27) with a Newton type algorithm implies to inverse a $(d+1)\times(d+1)$ matrix. In addition, for each $\theta$, obtaining the partial derivatives of (27) requires the computation of $(K_{x,x}^{\theta}+B)^{-1}$ [67]. Thus, at each step of the algorithm the complexity is $O(n^{\prime 3}+d^{3})$. Assuming the starting point $\theta_{\operatorname*{start}}$ is close to an optimal $\theta^{*}$, based on the same analysis as in section 4.1.3, the complexity to find $\theta^{*}$ is

[TABLE]

with $L$ the Hessian Lipschitz constant.

Finally, obtaining $\hat{q}_{\tau}^{\theta^{*}}$ from (29) implies inversing the matrix $K_{x,x}^{\theta}+B$ that is in $O(n^{\prime 3})$. So the whole complexity is

[TABLE]

5.2 Bayesian variational regression

Quantile kriging requires repeated observations to obtain direct observations of the quantile and make the hypothesis of Gaussian errors acceptable. Variational approaches allow us to remove this critical constraint, while setting a more realistic statistical hypothesis on $\varepsilon$. Starting from the decomposition of Eq.23, $\varepsilon(x)$ is now assumed to follow a Laplace asymmetric distribution [98, 55], implying:

[TABLE]

with the priors on $q$ and $\sigma$ that has to be fixed. Note that such assumption may be justified by the fact that at $\sigma$ fixed, minimizing the empirical risk associated to the pinball loss is equivalent to maximizing the asymmetric Laplace likelihood, which is given by

[TABLE]

According to the Bayes formula, the posterior on $q_{\tau}$ can be written as

[TABLE]

As the normalizing constant is independent of $q_{\tau}$, considering only the likelihood and the prior is enough. We obtain

[TABLE]

Because of the Laplace asymetric likelihood, contrary to the classical kriging model, here the posterior distribution (32) is not Gaussian anymore. Thus, it is not possible to provide an analytical expression for the regression model. To overcome this problem, [10] used a variational approach with an expectation-propagation (EP) algorithm [60], while [1] used a variational expectation maximization (EM) algorithm which was found to perform slightly better.

5.2.1 Variational EM algorithm

The EM algorithm was introduced in [26] to compute maximum-likelihood estimates. Since then, it has been widely used in a large variety of fields (see [57] for more details). Classically, the purpose of the EM algorithm is to find $\zeta$ a vector of parameters that define the model and that maximizes $p(\mathcal{Y}_{n}|\zeta)$ thanks to the introduction of the hidden variables $z=(z_{1},\cdots,z_{M})$. However dealing with this classical formalism implies that $p(z|\mathcal{Y}_{n},\zeta)$ is known or some sufficient statistics can be computed (see [90, 69] for more details), which is not always possible. Using the variational EM framework is a possibility to bypass this requirement [90], by approximating $p(z|\mathcal{Y}_{n},\zeta)$ by a probability distribution $\tilde{p}$ that factorizes under the form

[TABLE]

Starting from the log-likelihood $\log p(\mathcal{Y}_{n}|\zeta)$, thanks to Jensen’s inequality, it is possible to show that:

[TABLE]

where

[TABLE]

and $\operatorname*{kl}$ is the Kullback-Leibler divergence:

[TABLE]

As presented on Figure 4, the EM algorithm can be viewed as a two-step optimization technique. The lower bound $\mathcal{L}$ is first maximized with respect to $\tilde{p}$ (E-step) so that to minimize $\operatorname*{kl}(\tilde{p}||p)$. Next the likelihood is directly maximized with respect to the parameter $\zeta$ (M-step).

Classically and in the following, the E-step optimization problem is analytically solved (see [90] for details about the computation). In this particular case, the quantities $\tilde{p}_{j}$ are limited to conditionally conjugate exponential families. But note that different strategies have been developed to relax this assumption (generally it comes with an higher computational complexity) so that creating more flexible models (see the review [99] for instance).

5.2.2 Variational EM applied to quantile regression

Following [1], let us suppose that

[TABLE]

[TABLE]

with $\operatorname*{IG}$ defining the inverse gamma distribution and for sake of simplicity, $m=0$. Note that contrary to the formalism introduce with QK, here $\sigma$ is taken as a random variable. To allow analytical computation, let us introduce an alternative definition of the Laplace distribution [55, 46]:

[TABLE]

where $\mu_{y_{i}}=q_{\tau}(x_{i})+\frac{1-2\tau}{\tau(1-\tau)}\sigma w_{i}$, $\sigma_{y_{i}}=\frac{2}{\tau(1-\tau)}\sigma^{2}w_{i}$ and $w_{i}$ is distributed according to an exponential distribution of parameter $1$.

The distribution of $q_{\tau}$ at a new point $x_{*}$ is given by averaging the output of all Gaussian models with respect to the posterior $p(q_{\tau},\sigma,w|x,Y)$:

[TABLE]

Here the crux is to compute the posterior

[TABLE]

To do so, in [1] the authors use $z=(q_{\tau}(\mathcal{X}_{n}),w,\sigma)$ as hidden variables and $\zeta=\boldsymbol{\theta}\in\mathds{R}^{d+1}$ as parameters and the variational factorization approximation

[TABLE]

The EM algorithm provides a nice formalism here. Although the goal is to find $\boldsymbol{\theta}$ such that p\big{(}\mathcal{Y}_{n}|\mathcal{X}_{n},\boldsymbol{\theta}\big{)} is maximal, the algorithm estimates the underlying GP (i.e. p\big{(}q_{\tau},\sigma,w|\mathcal{D}_{n}\big{)}) that is able to have a likelihood as large as possible. Then the estimated value $p(q_{\tau},\sigma,w|\mathcal{D}_{n})$ is plugged into (34) to provide the final quantity of interest.

E-step.

Because the posteriors are conjugated, it is possible to obtain an analytical expression of the optimal distribution $\tilde{p}(q_{\tau})$:

[TABLE]

where

[TABLE]

and

[TABLE]

with

[TABLE]

The posterior on $w_{i}$ is a Generalized Inverse Gaussian GIG$(1/2,\alpha_{i},\beta_{i})$ with :

[TABLE]

and

[TABLE]

Due to numerical problems, in [1] the authors use the restriction $\tilde{p}(\sigma)=IG(a,b)$. Finding the best $a,b$ is done numerically. Then finding the best $a,b$ is equivalent to maximizing:

[TABLE]

with

[TABLE]

and

[TABLE]

M-step.

Ignoring terms that do not depend on $\boldsymbol{\theta}$, we obtain the lower bound:

[TABLE]

The optimization of $\mathcal{\tilde{L}}$ with respect to $\boldsymbol{\theta}$ is done using a numerical optimizer.

Recalling the goal is to compute (34), thanks to (35), we make the approximation:

[TABLE]

Then we obtain

[TABLE]

where

[TABLE]

Finally, as explained in section 5.1.1, the quantile estimator $\hat{q}_{\tau}$ is selected as $\bar{q}_{\tau}$. Algorithm 6 details the implementation of the VB method.

5.2.3 Computational complexity

E-step.

The complexity of this step is in

[TABLE]

In fact the algorithm computes $\Sigma_{\theta}$ that implies inversing a matrix of size $n\times n$.

M-step.

Optimizing $\mathcal{\tilde{L}}$ with a Newton type algorithm costs $O(n^{3}+d^{3})$ at each iteration (for details refer to the optimization description of (27)). Assuming the starting point $\theta_{\operatorname*{start}}$ is close to an optimal $\theta^{*}$, based on the same analysis as in section 4.1.3, the whole complexity is in

[TABLE]

Overall complexity.

At each iteration of the EM algorithm, the computation cost is $O(L(d^{3}+n^{3})+n^{3})$. The final complexity is obtained by multiplying by the number of iterations $n_{\operatorname*{it}}$ of the EM algorithm. Thus, the overall complexity is in

[TABLE]

6 Metamodel summary and implementation

In this section we detail our implementation procedure. After providing a summary of the six metamodels in Table 1, we present the packages we used and the hyperparameters we chose (which hyperparameters we set and which hyperparameters we optimized). We then describe the procedure we used to optimize the hyperparameters (optimization strategies and evaluation metrics).

6.1 Summary of the models

Table 1 lists the analytical expressions of the six metamodels, along with the associated underlying quantity.

6.2 Packages and hyperparameter choices

Each method depends on many parameters that can be tuned to improve performance, for example the choice of the kernel function and the value of its parameters for RK, QK and VB or the penalization factor for RK and NN. Here, to limit the computational burden, we chose to optimize only the most critical ones. When possible, for the other parameters, we applied the arbitrary choices and values made by the authors of the original papers. Most changes were made to improve robustness. Below, we describe our experimental settings, also summarized in Table 2.

Nearest Neighbors.

We set $d()$ as the Euclidean distance and optimized only the size $K$ of the neighborhood.

Random forest.

In this case, the only hyperparameter we optimized was the maximum size of the leaves $m_{s}\in\mathds{N}^{*}$. Regarding the number of trees, we noticed that the metamodel needs many more trees than are needed for the estimation of the expectation. In some problems, the metamodel needs up to 5,000 trees to stabilize the variance. Thus, in our experiments we set the number of trees at 10,000 in all cases. We set the number of dimensions considered for the split at the default choice $d/3$. The depth of the tree is not constrained and the splitting rule is based on Eq. 7. We used the R package QuantregForest [59].

Neural network.

Based on [18], we set the number of hidden layers at one and the transfer function as the sigmoid. The optimization algorithm is a Newton method, we set $E_{K}$ at $1/2^{K}$ with $K=1,2,5,10,15,20,25,30,35$ and the number of multistarts to optimize the empirical risk at five. We optimized the number of neurons $J_{1}$ in the hidden unit and the regularization parameter $\lambda\in\mathds{R_{+}}$. The metamodel is generated using the R package qrnn [18].

Regression in RKHS.

The kernel was set as a Matérn 5/2. We optimized the length scale parameters $\theta\in\mathds{R}_{+}^{d}$ and the regularization hyperparameter $\lambda\in\mathds{R_{+}}$. Optimization (20) is done with the quadratic optimizer quadprog [89].

Quantile Kriging.

The kernel was set as a Matérn 5/2. The number of repetitions was set according to the total number of observations (see Table 4). We optimized the length scale hyperparameter $\theta\in\mathds{R}_{+}^{d}$ and variance hyperparameter $\rho\in$ $\mathds{R}_{+}$. QK is implemented in the R package DiceKriging [72].

Variational regression.

The kernel was set as a Matérn 5/2, the number of EM iterations $n_{\operatorname*{it}}$ at $50$. We optimized the length scale hyperparameter $\theta\in\mathds{R}_{+}^{d}$ and variance hyperparameter $\rho\in$ $\mathds{R}_{+}$. The implementation is based on the Matlab code provided in [1].

6.3 Tuning the hyperparameters

In the previous section, we defined the hyperparameters we wanted to optimize for each method. In fact, once the type of metamodel is chosen, the quantile estimator is given by a function $\hat{q}_{\Theta}:\mathcal{X}\rightarrow\mathds{R}$ where $\Theta\in\mathds{R}^{v}$ are called hyperparameters and $v$ is metamodel dependent. Hyperparameter optimization (also known as model selection) is an essential procedure when dealing with non-parametric estimators. Although $\hat{q}_{\Theta}$ may be very efficient on $\mathcal{D}_{n}$, the prediction may perform very poorly on an independent dataset $\mathcal{D}^{\prime}_{p}$. The goal is to find the ${\Theta}$ that provides the best possible generalized estimator. In the following, we present the validation metric used to optimize the hyperparameter values and detail the hyperparameters optimization procedure associated with each method.

6.3.1 Metrics

In the standard conditional expectation estimation, the validation and performance metrics are both based on $\left\|\hat{m}_{\Theta}-m\right\|_{L^{2}}$, where $\hat{m}_{\Theta}$ is the estimator and $m$ the targeted value. With the quantile estimation procedure the two metrics are no longer the same. The goal is to find $\hat{q}_{\Theta}$ such that $\left\|\hat{q}_{\Theta}-q\right\|_{L^{2}}$ is as small as possible. However, $q$ is unobserved so the validation metric cannot be based on the $L^{2}$ norm. Here we present two metrics able to measure the generalization capacity of a quantile metamodel.

Bayesian metamodels (QK and VB) have their own validation metric, this is the likelihood function that can be maximized with respect to $\Theta$. For quantile kriging, we use (27) while in the variational approach we use (36). The optimal hyperparameters are then:

[TABLE]

The second metric available for all metamodels is $k$-fold cross-validation associated with the pinball loss. The metric can be computed as follows. First, the data are split into $k$ parts, then the model is trained on $\mathcal{D}_{-j}$, the training set without the $j$-th part. The performance is evaluated on the remaining part $\mathcal{D}_{j}$. As the quantile minimizes the pinball loss (on $\mathcal{D}_{j}$), the evaluation metric is

[TABLE]

where $n_{j}^{\prime}$ is the number of observations in each fold. The optimal cross-validation hyperparameters are then:

[TABLE]

In our experiments, we chose $k=5$ to limit the computational cost. However, we observed empirically that choosing a larger $k$ did not substantially modify the performances of the metamodels. Although cross-validation is available for QK, we chose to stay in the spirit of the methods and to only use likelihood to select hyperparameters. Our choice is supported by [5] that does not show a clear improvement using cross validation instead of maximum likelihood techniques.

6.3.2 Hyperparameters optimization procedure

Both likelihood functions come with analytical derivatives, enabling the use of gradient-based algorithms. However, since both functions are multi-modal, multi-start techniques are necessary (and generally very efficient, see [38]) to avoid being trapped in local optima. To account for the increasing difficulty of the optimization task with the dimension while limiting the computational cost, we ran $n_{\operatorname*{start}}=20d$ optimization procedures from different starting points $\theta_{\operatorname*{start}}$ and chose the set of starting points based on a maximin Latin hypercube design.

For QK, the BFGS algorithm is used to optimize (27). For VB, two derivative-based optimizers are used alternately for the E- and M-steps. Since each step may lead the algorithm toward a local minimum, we chose to apply the multi-start strategy in the entire EM procedure.

Optimization of the cross-validation metric (38) is done under the black box framework, since no structural, derivative or even regularity information is available. Hence, all optimizations are carried out using the branch-and-bound algorithm named Simultaneous Optimistic Optimization (SOO) [62]. SOO is a global optimizer, hence robust to local minima.

We used [87] to parallelize the computations.

6.3.3 Oracle metamodels

Each method presented in this paper is a trade-off between power and the difficulty of finding good hyperparameters. A good method should be powerful (i.e. provide flexible fits) but easy to tune. In order to assess the strengths and weaknesses of the hyperparameter tuning methods in addition to standard metamodels, we provide what we call the oracle metamodels for each problem. Instead of using the cross-validation or likelihood metric, the oracle tunings are directly based on the evaluation metric

[TABLE]

where $n_{\text{test}}$ is the size of the test set. In a sense, they provide a upper bound on the performance of each method. This allows us to show which metamodels have the potential to tackle the problems and which are intrinsically too rigid or make poor use of information. In addition, this allows us to directly assess the quality of the validation procedure.

7 Benchmark design and experimental setting

Many factors can affect the efficiency of methods to estimate the right quantile. For our benchmark system, we considered five models or test cases to evaluate the performance of the six metamodels. We decided to focus primarily on the dimensionality of the problem, the number of training points available, the signal-to-noise ratio defined as

[TABLE]

and the pdf value at the targeted quantile for test cases in which the distribution shape and the distribution spread (i.e. level of heteroscedasticity) can vary significantly. Our two objectives were to:

discover if there is a single best method for all factors variations considered or specific choices depending on the configuration at hand, and 2. 2.

assess the performance of the quantile regression, and in particular, the configurations for which the current state-of-the-art is insufficient.

A full 3D factorial experimental design was used to analyze the efficiency of the metamodels, the 3 factors being the test case (4 test cases), the number of training points (4 levels) and the quantile order (0.1, 0.5 and 0.9). We used part of this complete design to focus our analyses on the characteristics of the test cases (dimension, pdf shape and heteroscedasticity).

7.1 Test cases and numerical experiments

Test case 1

is a $1$d toy problem on $\left[-1,1\right]$ defined as

[TABLE]

with $\xi=\eta\mathbb{1}_{\eta\leq 0}+7\eta\mathbb{1}_{\eta>0}\text{ where }\eta\sim\mathcal{N}(0,1).$

The signal-to-noise ratio is $\operatorname*{SNR}\approx 0.5$, it is consider as small. The pdf value $f(x,q_{\tau})$ varies substantially according to $x$ for all $q_{\tau}$. Indeed, on the interval $[-0.5,-0.3]$, for all values of $q_{\tau}$, the pdf is very large (almost equal to $+\infty$) because the variance of the distribution on this interval is very small. In contrast, for $\tau=0.9$ (resp. $\tau=0.1$) the pdf is very small in the interval $[0.6,1]$ (resp. $[-1,-0.6]$). In $[0.6,1]$ the pdf of the $0.9$-quantile is equal to the pdf of the $0.9$-quantile of a normal distribution of variance $49(0.2+3x^{3})^{2}$ that is, for example, approximately equal to $0.01$ for $x=0.9$. Because of this important variations, we consider the values of the pdf for all quantiles of interest as variable.

Test case 2

is a $2$d toy problem on $\left[-5,5\right]\times\left[-3,3\right]$ based on the Griewank function [28], defined as

[TABLE]

with

[TABLE]

and $\xi=\eta\mathbb{1}_{\eta\leq 0}+5\eta\mathbb{1}_{\eta>0}\text{ where }\eta\sim\mathcal{N}(0,1).$

The signal-to-noise ratio is $\operatorname*{SNR}=0$, it is consider as small. The pdf value $f(x,q_{\tau})$ varies substantially according to $q_{\tau}$. Indeed, for $\tau=0.1$ (resp. $\tau=0.9$) the pdf is small (resp. very small), more precisely the pdf of the $0.1$-quantile (resp. $0.9$-quantile) is equivalent to the pdf of the $0.1$-quantile (resp. $0.9$-quantile) of a normal distribution of variance $G^{2}$ (resp. $25G^{2}$) with $G$ that takes values in $[0,2]$. That implies $f(x,q_{\tau})$ varies with respect to $x$ as well. Note that close to $x=(0,0)$ the pdf is very large because of the very small variance of the associated distribution. The pdf at the $0.5$-quantile is discontinuous. To the left of the $0.5$-quantile the pdf is equal to the pdf at the median of a normal distribution of variance $G^{2}$ while to the right it is equal to the pdf at the median of a normal distribution of variance $25G^{2}$. To treat this discontinuity, the value of the pdf is considered as the maximum of both sides. Thus, for this particular case the pdf is considered as large.

Test case 3

is a $1$d toy problem based on the Michalewicz function [28] on $\left[0,4\right]$, defined as

[TABLE]

with $\xi=3\eta\mathbb{1}_{\eta\leq 0}+6\eta\mathbb{1}_{\eta>0}\text{ where }\eta\sim\mathcal{N}(0,1)$.

The signal-to-noise ratio is $\operatorname*{SNR}\approx 0.04$, we consider it as small. The pdf value $f(x,q_{\tau})$ varies substantially according to $q_{\tau}$ and $x$. The conditional distribution of this problem is not a classical one but it is close to the distribution of $-\mathcal{X}^{2}$ with one degree of freedom. It implies at $x$ fixed, the pdf value increases with $\tau$. For the $0.1$-quantile (resp. $0.9$-quantile) the mean of the pdf in the interval $[0,2.5]$ is approximately $0.05$ (resp. $9.1$). According to $x$ the pdf varies as well. For $x\in[3.5,4]$ the variance of the conditional distribution is very small thus the pdf near the quantiles of interest is very large. The value of the pdf at the $0.5$-quantile is between the value of the pdf at the $0.1$ and $0.9$ quantile. Thus, we consider the pdf at $q_{0.1}$ as globally small and at $q_{0.5}$ and at $q_{0.9}$ as globally large.

Test case 4

is a $9$d toy problem based on the Ackley function [3] on $\left[-1,-0.7\right]\times\left[0,1\right]\times\left[-0.7,-0.3\right]\times\left[0.5,1\right]\times\left[-1,-0.5\right]\times\left[-3,-2.6\right]\times\left[-0.1,0\right]\times\left[0,0.1\right]\times\left[0,0.8\right]$, defined as a function

[TABLE]

with

[TABLE]

and

[TABLE]

with $a=10$, $b=2\times 10^{-4}$, $c=0.9\pi$ and $\xi$ follows a log-normal distribution of parameters $(0,1)$.

The signal-to-noise ratio is $\operatorname*{SNR}\approx 2.3$, it is consider as large. The pdf value $f(x,q_{\tau})$ varies substantially according to $q_{\tau}$ and $x$. Indeed, for $\tau=0.1$ (resp. $\tau=0.9$) the pdf is large (resp. very small). At the conditional distribution of this test case is a log-normal distribution of parameters \big{(}\log R(x),1\big{)}. The expectations of the pdf value at different values of $\tau$ are $\mathds{E}(f(.,q_{0.1}))=0.1$, $\mathds{E}(f(.,q_{0.5}))=0.04$ and $\mathds{E}(f(.,q_{0.9}))=0.005$. Thus, the pdf at $q_{0.1}$ is consider as large, the pdf at $q_{0.5}$ and $q_{0.9}$ are consider as small and very small.

To provide a better intuition about this problem, we plotted what we call the marginals. For all dimensions except the $j$-th, the values of the input are fixed to $x_{-j}\in\mathds{R}^{8}$ and the $j$-th dimension varies. Figure 7 represents the evolution of the quantiles w.r.t. the $j$-th dimension for two different $x_{-j}$ and for $j=1,\cdots,9$. In particular it shows that the difference between the $0.1$ and $0.9$-quantile depends significantly on $x$.

Note that on those four toy problems, the random term $\xi$ is defined such that the resulting distribution of $Y$ would be strongly asymmetric. As $\xi$ is also multiplied by a factor that depends on $x,$ the distribution of $Y$ is also heteroscedastic. The first three toy problems are represented in Figure 5 and some illustrations of the fourth test case are available Figure 7.

Test case 5

is based on the agronomical model SUNFLO, a process-based model which was developed to simulate sunflower grain yield (in tons per hectare) as a function of climatic time series, environment (soil and climate), management practices and genetic diversity. The full description of the model is available in [19]. In the regression model we consider $\mathcal{X}$ corresponding to nine macroscopic traits that characterize the sunflower variety. Although SUNFLO is a deterministic model, for each simulation the climatic time series are randomly chosen within a database containing $190$ years of weather records, which makes the output stochastic (see also [63] for more details).

The signal-to-noise ratio is $\operatorname*{SNR}\approx 0.1$. The pdf value $f(x,q_{\tau})$ varies substantially according to $q_{\tau}$, more precisely $\mathds{E}(f(.,q_{0.1}))=0.17$, $\mathds{E}(f(.,q_{0.5}))=0.15$ and $\mathds{E}(f(.,q_{0.9}))=0.05$. Thus, the pdf at $q_{0.1}$ and at $q_{0.5}$ are consider as large, and the pdf at $q_{0.9}$ is consider as small.

In addition, the shape of the distribution varies significantly with respect to the input as illustrated Figure 6.

Numerical experiments.

On all problems, we consider four sample sizes. Those sizes depend on the dimension and are empirically chosen so that the smallest size corresponds to the minimal information required by the metamodels to work and the largest size is chosen keeping in mind the potentially high cost of simulators. Besides, our focus is on computer experiments, where data sizes rarely exceed thousands of points. For the $1d$ problems, the points are generated on a uniform grid. For the $2d$ and $9d$ problems, the observations are taken on a Latin hypercube design optimized for a maximin criterion to ensure space-filling [32]. The same samples are used by all methods except QK, as it requires repetitions. For QK, the number of distinct points and number of repetitions depends on the budget. The different sample sizes are reported in Table 4. Finally, for each sample size and problem, 10 samples are drawn in order to assess robustness with respect to the data.

7.2 Structuration between the questions and the numerical setting

Factors.

Three factors are explicit in our benchmark system: the number of training points, the problem dimension and the quantile level. The other factors depend on the characteristics of the problem concerned: shape variation, pdf value at the quantile, level of heteroscedasticity, signal-to-noise ratio. For all four test cases, we consider three quantile levels: $0.1$, $0.5$ and $0.9$. Note that due to the asymmetry of the problems, learning for the $0.1$ and $0.9$ quantiles is not equivalent in terms of difficulty. Indeed, when the response is heteroscedastic (a variance/spread depending on $x$) and/or when the shape of $\mathds{P}_{x}$ varies in $x$, the pdf $f(x,q_{\tau})$ may vary in $x$ as well. Intuitively, quantiles with large pdf values are easier to learn, as the data points may be closer to them. Figure 8 illustrates this effect. Table 3 summarizes the characteristics of our design concerning the number of training points with respect to the dimension of the test case. To make our results easier to analyze, we divided the problems into groups that allow us to focus on subsets of factors.

Focus 1: is there a universal winner?

To provide a universal ranking of the methods, we use all test cases, training points and quantile levels. As highlighted in Table 3, we created a set of different problems representative of a large number of characteristics that could be met dealing with any quantile regression problem. Note that our benchmark system is slightly biased towards small-dimensional problems, since only three-fifths of the cases have a dimension higher than two.

Focus 2: what behavior with respect to the dimension, the number of training points, signal-to-noise ratio and pdf value?

To analyze the effects of these factors on the performance of the methods, we combine toy problems $1$, $2$, $3$, $4$ and the SUNFLO model. Note that once the pdf value is taking as a explanation variable, toy problem $1$ is excluded from the group because the pdf value near all the studied quantiles cannot be classified as large or small.

7.3 Performance evaluation and comparison metrics

Assessing the performance of quantile regression is not an easy task when only limited data are available. Here, since we are considering toy problems (exept for SUNFLO, in that case the true quantile values are taken as the quantiles of the $190$ years of weather records), the true quantile values can be approximated with precision, so we can evaluate the $L^{2}$ error for each emulator. The value of the criterion is privided by (39).

We chose $n_{\text{test}}=250$ for the $1d$ problems and $n_{\text{test}}=4000$ for the others. Now, since the problems vary in difficulty and in their response range (Figure 5), $E_{L^{2}}(\hat{q})$ cannot be aggregated directly over several problems or configurations. To do so, we normalize this error by the error obtained by a constant model (the constant being taken as the quantile of the sample):

[TABLE]

where $\operatorname*{CQ}$ stands for constant quantile.

As an alternative criterion, we consider the ranks of the metamodels based on their $L^{2}$ error. Although ranks do not provide information regarding the range of errors, they are insensitive to scaling issues, which makes aggregation between configurations more sensible. They allow us to assess whereas any method consistently outperforms others, regardless of overall performance.

8 Results

8.1 Focus 1: overall performance and ranks

First, we consider the overall performance and ranks, integrated over all runs. We have considered $5$ test cases, for each test case we have generated $4$ sizes of training sets, for each test case and sample size 10 samples are drawn and for each of this occurrences we have estimated $3$ different conditional quantiles. Thus, we have $5\times 4\times 10\times 3=600$ experiments. They are shown as boxplots in Figure 9. Based on these ranks (Figure 9, left), VB appears to be the best solution since it is ranked either 1st or 2nd in 50% of the problems. RK is in second position, its median is the same as VB but it is generally ranked between 2nd an 3rd. In addition RK seems slightly less risky than VB in the sens that it is almost never rank 5th or 6th. KN is the worst since its median rank is equal to five.

However, all boxplots range between 1 and 6, indicating that no method is outperformed by another on all problems. This finding is reinforced by the performance boxplots (Figure 9, right), where all median performances are similar (VB and RK being the best and QK the worst which means QK may be very bad sometimes), and the variance is very large. Indeed, the errors range from 2% for NN (of the error achieved by a constant metamodel) to 500% for QK, all methods experiencing cases with more than 100% error (i.e. situations where they are worse than the constant metamodel).

8.2 Focus 2: dimension, number of training points and pdf value

8.2.1 Performance according to the constant quantile

In this section, we analyze the performance of the methods with respect to the pdf value and the number of points.

Sample size:

Figure 10 shows the performances of the methods grouped according to the size of the sample. As expected, the performances increase with the size of the sample. For size $1$ ($n\approx 50d$), the distribution of $E_{\operatorname*{CQ}}$ of all the metamodels is almost centered around 100%, implying that these correspond to limit cases for quantile regression since the metamodels do not outperform the constant metamodel (although in some cases the error is as small as 40%). For size $4$ ($n\approx 300d$), the median performance is roughly 50% (twice as accurate as the constant metamodel). BV, RK and especially NN experience situations with very accurate models. However, all the methods also experience bad performances (error greater than 100%) in the large sample regime. Unfortunately, from Figure 10 we can conclude that no method is sufficiently robust in all cases.

Signal-to-noise ratio.

Figure 12 groups performance with respect to the signal-to-noise ratio. According to the figure the performance depends to a great extent on the signal-to-noise ratio. Dealing with a high signal and considering the performance, there is a clear difference between the statistical methods (KN, RF) and the four others (RK, NN, VB, QK) while this difference is not visible in the small signal setting. The impact of the signal-to-noise ratio is strong on the performance. Dealing with a high signal, the median of the performance for NN, RK, QK, VB is close to 20% with a small variance while if the signal is small, the error is larger (the median is above 50% for all methods) and the variance is larger.

In the following we focus our interest only on cases with small signal-to-noise ratio. Indeed as the performance of RK, NN, QK and VB are close to each others, we think we do not have enough experiments to extract patterns.

Pdf value with small signal-to-noise ratio.

Figure 12 groups performance with respect to sample size (either small, i.e. level 1 and 2 or large, i.e. level 3 and 4) and pdf value (according to Table 3). According to the Figure 11, the performance depends to a great extent on the pdf value in the neighborhood of the targeted quantile. With a small $n$, the pdf value has no significant impact on the median of the performance (except for VB) but it does have an impact on the lower bound of the error. More precisely, the median of the error does not depend on the pdf value in the case of a small pdf but sometimes the metamodel errors are sensibly smaller when the pdf is large. With large samples, both the median and the lower bound of the error depend on the pdf value. Metamodels may be very good when the pdf is large, for example $20$ times better than the constant metamodel for NN whereas the error appears to have a lower bound when the pdf is small even with large $n$. In addition, for a problem with a small $n$ and a large pdf, the performance is similar to the performance for problems with a large $n$ and a small pdf (Figure 12, see the two columns in the center).

8.2.2 Rank in the context of a low signal-to-noise ratio

Pdf value.

Figure 13 shows clearly that when the pdf is large, VB is the best model while when the pdf is small (and the problem is heteroscedastic), VB is less good than RN, RK and KN. This observation is supported by Figure 12 which reveals a strong contrast between the performance of the VB method. QK is poor in both cases, whereas RK performs comparatively better with small pdf.

Sample size.

Figure 14 shows that the number of points has a major impact on the ranking of some methods. The ranking of QK and VB is relatively insensitive to the size of the sample. The other methods are less distinguishable when the sample size is small than when it is large. With a small sample KN, RF, NN and RK are comparable, whereas when the sample size increases, NN and RK clearly outperform KN and RF. For the largest size, NN is slightly better than VB.

Dimension.

Figure 15 groups performance based on dimension. The first contrast is the permutation between RK and NN. With a small dimension, RK is better than NN but the relative performance of NN increases w.r.t. the dimension. With small dimensions, RF and KN are comparable, but with high dimensions, RF outperforms KN.

High dimension, small pdf.

Figure 16 shows an extreme case in which the pdf is small but the dimension is high. With a small $n$, the best method is clearly RF followed by RK and KN. VB and QK are not well ranked. With a larger $n$, as mentioned above, NN and VB are better but with large variance, while RF, RK and KN rank less well.

9 Extensions and open questions

9.1 Effect of hyperparameter tuning

In the following we define

[TABLE]

the performance gap between the regular metamodel and its oracle performance (the loss in performance between actual hyperparameter tuning and oracle tuning). Figure 17 gives the average values of $\Delta E$ aggregated respectively over all problems and only aggregated over the problems with a large pdf, and considering the effect of dimension and sample size. In high dimension, the easiest methods to tune are KN, NN, and RF. In our study, KN and RF have a single hyperparameter to tune regardless of the dimension, and NN has two. This is clearly an advantage in terms of robustness in high dimension. The other methods are kernel-based and require the tuning of at least $d+1$ hyperparameters. This consistently affects RK and QK, but affects VB only in the case of small pdf, while it is the most stable method in the other cases.

With a small dimension, all the methods have roughly the same number of hyperparameters. The most noticeable change compared to the case of a high dimension is the good performance of RK, while NN becomes comparatively the most difficult method to train.

9.2 On the methods’ behavior

Statistical order methods.

As presented on Figure 11 this methods are not relevant on high signal-to-noise regime. A possible explanation is the difficulty to fit smooth variations with high amplitudes with a piecewise constant model.

Concerning the low signal-to-noise regime, it is clear from Figure 15 that KN performs poorly in a high dimension. This may be due to the irrelevance of the Euclidean distance when there is a significant increase in dimension. RF clearly outperforms KN in this situation, as it is able to produce better neighborhoods than the Euclidean distance. Overall, (compared with the other methods), RF performance increases with dimension. This may be due to the fact that it has fewer hyperparameters to tune.

Functional methods.

As presented on Figure 11 this methods are relevant on high signal-to-noise regime.

Concerning the low signal-to-noise regime, Figure 14 shows that NN works poorly in a small sample setting, but it is one of the best methods when the sample is large. This result reflects the high flexibility of NN. Too much flexibility leads to overfitting when the sample is small. In contrast, when the number of points is large, NNs are able to fit the data very well (e.g. Figure 10, Size $4$). According to Figures 10 and 14, RK is a robust method. Its robustness in both small and large data settings can be attributed in part to the selected kernel. If the selected kernel is sufficiently smooth (here continuous and derivable), the resulting metamodel cannot produce instable results. However, it seems (Figure 10, Size $3$ and $4$) that this lack of flexibility may affect the performance with an increase in the size of the data set. In this case, more flexible methods like NN may outperform RK. The contrast between RK and NN shown in Figure 15 can be explained by the level of difficulty associated with each method involved in finding good hyperparameters (as explained above).

Bayesian models.

As presented on Figure 11 this methods are relevant on high signal-to-noise regime.

Dealing with low signal-to-noise ratio the QK method under-performs comparing to others. One possible explanation is the erroneous assumption in Equation (28) that the noise is centered, which is more critical for extreme quantiles. Another possible explanation lies in the small number of replica. The local inference (that uses statistical orders) is biased and in a low signal-to-noise regime it has high variance. In addition, the increasingly bad performance of QK with an increase in dimension (Figure 10) is likely a consequence of the fact that empty areas become larger in high dimensions.

VB is one of the best methods presented in our paper. Figures 13 and 17 show that VB is also the most dependent on the pdf value. When the pdf is large, it is the best method whereas when the pdf is small (and the shape of the distribution and/or the variance depend on the input), it may be the worst. The explanation lies in the philosophy of the model. In the case of RN and RK, the model complexity (i.e. smoothness) is almost entirely related to the regularity of parameter $\lambda$ that is selected by cross-validation. Hence, the model cannot excessively overfit and cannot perform very poorly. With Bayesian methods, the regularization is included in the model hypothesis: in our setting, the quantile is assumed to be a Gaussian process with covariance function $k_{\theta}(.,.)$, so $\theta$ performs the regularization. We observed that if the local quantity of information (roughly the product of the number of points times the pdf value in the neighborhood of the quantile) is too small comparing to the information needed to fit well the quantile, the metamodel tends to interpolate the available data. When sufficient information is available, the optimization of the marginal likelihood provides $\theta$ values that allow a good trade-off between flexibility and smoothness. This is likely the reason why VB is easily beaten by RN and RK when the pdf is small.

9.3 Varying shape and heteroscedasticity.

If the shape of $\mathds{P}_{x}(Y)$ or the variance of $Y_{x}$ (heteroscedasticity) vary w.r.t. $x$, then $f(x,q_{\tau})$ may vary in $x$. Figure 18 illustrates the ability of RF, RK and VB to estimate quantiles of a distribution with a strongly varying shape. In this problem, as depicted in Figure 18 (top row), the quantiles are not perfectly estimated but the metamodels provide good indications about the shape of the true distribution. However, as can be seen in Figure 18 (bottom row), the methods can present strong instabilities. Here, for a sample virtually indistinguishable from the one leading to accurate estimates, the median estimates largely overestimate the true values for large $x$ values. Such instabilities can be partly imputed to the difficulty of the task. However, this is also because no method is actually designed to deal with strongly varying pdf, as we explain below.

An ideal method would be almost interpolant for a very large pdf but only loosely fit the data when the pdf is small. Indeed, if the pdf is very large then the output is almost deterministic, thus the metamodel should be as closer as possible to the data. In the small pdf case, a point does not provide a lot of information. Information should be extracted from a group of points, that means the metamodel must not interpolate the data. However, most of the methods presented here rely on a single hyperparameter to tune the trade-off between data fitting and generalization: the number of neighbors for KN, the maximum size of the leafs for RF and the penalization factor for NN and RK. As a result, the selected hyperparameters are the ones that are best on average. Theoretically, this is not the case for the Bayesian approaches: QK accounts for it via the error variance $\sigma_{i}^{2}$ computed by bootstrap, and the weights $w_{i}$ (Eq. 33) allow VB to attribute different ”confidence levels” to the observations. However in practice, both methods fail to tune the values accurately, as we illustrate below. Figure 19 shows the three quantiles of toy problem 3 and their corresponding RF, RK and VB estimates. For $\tau=0.1$ in particular, the pdf ranges from very small ($x$ close to 0) to very large ($x$ close to 4). Here, RF and RK use a trade-off that globally captures the trend of the quantile, but cannot capture the small hill in the case of large $x$. Inversely, VB perfectly fits this region but dramatic overfitting occurs on the rest of the domain.

Finally, Figure 20 illustrates that this is not an issue of hyperparameter tuning. For each method, we show the oracle estimate, a tuning that tends to underfit and another that tends to overfit. One can see that no tuning is entirely satisfactory, since capturing the region with high pdf leads to overfitting on the rest of the domain and vice-versa.

We believe that further research is necessary to obtain estimators that intrinsically account for strong heteroscedasticity and varying shape. One possible direction is the use of stacking, in the spirit of [78]. Under the stacking framework the final estimator could be

[TABLE]

where $\{\hat{q}_{\theta_{i}}\}_{1\leq i\leq N}$ is a set of metamodel and $\{g_{i}(x)\}_{1\leq i\leq N}$ is a set of weight functions. Choosing $\{\hat{q}_{\theta_{i}}\}_{1\leq i\leq N}$ such that they correspond to different pdf values might provide more flexible estimates.

9.4 On the non-crossing of the quantile functions

While the quantile functions (for different quantile levels) may obviously never cross, unfortunately, their estimators may not always satisfy this property. This is a well-known issue against which none of the methods presented here is immune.

The neighborhood approaches first estimate the CDF, then extract the quantiles. If the hyperparameters are the same for all quantiles, crossing is impossible. However in our setting, different neighborhood sizes were used for different quantiles.

With functional analysis approaches, crossing may happen even if each quantile is built with the same hyperparameters. In the literature, authors have produced methods to address this issue. It could be reduced by the introduction of additional constraints in the model [86] or by the construction of a new model that intrinsically produces non-crossing curves [73]. However in both cases, the dimension of the optimization problem then increases significantly.

Finally the stochastic process approaches estimate each quantile in independent Gaussian processes, so crossing may occur. While the number of training points is larger than what we consider here ($10^{4}$ repetitions for each input point), [17] use GP and takes into account all the quantile orders at the same time and thus ensures non-crossing.

Another approach available for all the methods presented here is related to the rearrangement of curves or isotonic regression [6, 2]. The idea is to perform many quantile regressions with a large number of different values of $\tau$ or a large set of bootstrapped versions of $\mathcal{D}_{n}$ and then to rearrange the curves, thereby obtaining the whole distribution and then extracting the quantiles that by definition do not cross.

In theory, adding non-crossing constraints and predicting several quantiles simultaneously could improve the quality of the estimates (in particular as it might add some robustness). However, in practice, it also makes the model more rigid (i.e. a single regularization hyperparameter for all quantiles), and preliminary experiments have shown no gain in accuracy compared to independent predictions, despite a considerably higher computational cost. Hence, multi-quantile predictors were not considered in our study.

9.5 Assessment of prediction accuracy

To assess the accuracy of the results and the confidence that we can have in the estimation it could be useful to provide confidence intervals for the predictor. From this point of view the methods are not equal. With Bayesian approaches, theoretical confidence intervals are provided with the models. More precisely as the output model is Gaussian and as it returns the mean and the variance, confidence intervals can be created. For example the $0.9$-confidence interval is provided by

[TABLE]

However as presented on Figure 21, while the confidence intervals obtained from QK seem useful, the VB model is clearly overconfident.

The statistical order methods consider that inside each neighborhood the samples are i.i.d. Based on that, it is possible to use Wilks’ formula (see [68] for instance) or deviation inequality as presented in [88] to extract confidence intervals. For example, using Chernoff’s inequality, for any $\eta>0$, the confidence interval of order $1-\eta$ is as $\operatorname*{CI}(x)=[L_{K}(x),U_{K}(x)]$ with

[TABLE]

and

[TABLE]

Using this method enables the confidence intervals to be data dependent. For example in Figure 21 the confidence intervals are not symmetric. Note that while the confidence intervals obtained with this technique come with theoretical guarantees, they are very conservative and they depend a lot on the size of the training set.

Finally for all methods it is possible to use a bootstrap technique to create different regression models and then to aggregate them in order to create empirical confidence intervals. For example Figure 21 shows confidence intervals using a bootstrap technique for RK. The main drawback of such method is its computational cost.

To the best of our knowledge there are no other methods that sensibly improve this results. Thus, there is still a room for improvement concerning quantile regression model assessment.

10 Summary and perspectives

10.1 General recommendations

In this presentation we have introduced six metamodels for quantile regression. In the first part of the paper we have provided a full description of the six metamodels, first focusing on their theoretical basis, then discussing their implementation procedure. This part of the paper have enabled us to highlight the similarities and differences of the methods so that providing critical perspectives on the state of the art. The second part of the paper focused on performance comparison according to the dimension of the problem, the size of the learning set, the signal-to-noise ratio and the value of the pdf at the targeted quantile. We have compared the presented methods on $4$ toy problems in dimension 1, 2, 9 and on an agronomic model in dimension 9.

Figure 22 summarizes our findings. In a nutshell, when the signal-to-noise ratio is high RK, VB, QK and NN shows good results in our experimental setting but as soon as the signal-to-noise ratio decreases, quantile regression requires larger budgets and comparing the methods seems to be more complicated. Indeed, while the rule-of-thumb for computer experiments is a budget (i.e. number of experiments) $10$ times the dimension (see [52] for instance), as we work on problems with low signal-to-noise ratio, we found that no method was able to provide a relevant quantile estimate with a number of observations less than 50 times the dimension. For larger budgets, no method works uniformly better than any other. NN and VB are best when the budget is large. When the budget is smaller, RF, RK, KN are best when the pdf is small in the neighborhood of the quantile, in other words, when little information is available. However, VB outperforms all the other methods when more information is available, that is, when the pdf is large in the neighborhood of the quantile.

10.2 Possible ways of improvement

In our benchmark, we generally followed the approaches as presented by their authors. However, most of them could be improved. The optimization scheme of NN is the computational bottleneck of the method, which makes it the most expensive method in our benchmark system. One possible improvement would be using the BFGS algorithm (see [49] for details about the BFGS algorithm applied to non-smooth functions) or the ADAM algorithm [44] to optimize directly the empirical risk associated to the pinball loss. A faster scheme would allow more restarts, and hence improve robustness.

Another improvement concerns the splitting criterion (7) of RF, which is not designed for the quantile but for the expectation. This could lead to poor estimates for problems where quantiles are weakly correlated with expectations. Defining an appropriate splitting criterion could significantly improve the performance of this method.

In our experiments, QK used a predefined number of sampling points that were heuristically defined as a trade-off between space-filling and pointwise quantile estimation accuracy. The performance of QK could be significantly improved by optimally tuning the ratio between the number of points and repetitions, in the spirit of [8].

The KN method can naturally be extended to a variant that uses the whole sample instead of the $K$ nearest points. The weights associated to each point of the sample could be based on Gaussian or triangular kernel for example. This idea has been developed in [96] to estimate the conditional expectation but we think that it could be possible to adapt this approach to the estimation of the conditional quantile.

Finally, in practice, finding the best hyperparameters was the most difficult part of the proposed benchmark system. While this aspect is often toned down by authors, we believe hyperparameters tuning is a key practical aspect that remains a challenging problem in quantile regression.

Acknowledgments

This work is part of a Ph.D. of L. Torossian funded by INRA and Région Occitanie. The authors would like to thank Edouard Pauwels for the discussions on the topic of convex and non-convex optimization and Pierre Casadebaig for its suggestions and advices concerning the crop model sunflo.