# Adaptive finite-dimensional approximation of constrained Gaussian processes for large datasets

Hassan Maatouk, Didier Rullière, Xavier Bay

<figure><img src="image_1.png" /></figure>

## To cite this version:

Hassan Maatouk, Didier Rullière, Xavier Bay. Adaptive finite-dimensional approximation of constrained Gaussian processes for large datasets. 2025. hal-05175025

HAL Id: hal-05175025

https://hal.science/hal-05175025v1

Preprint submitted on 4 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Adaptive finite-dimensional approximation of
constrained Gaussian processes for large datasets

Hassan Maatouk¹*, Didier Rullière² and Xavier Bay²

¹LAMPS, Université de Perpignan Via Domitia, 52 Av. Paul Alduy, Perpignan,
66860 Cedex 9, France.

²Mines Saint-Étienne, Univ Clermont Auvergne, CNRS, UMR 6158 LIMOS,
Institut Henri Fayol, Saint-Étienne, F-42023, France.

*Corresponding author(s). E-mail(s): hassan.maatouk@univ-perp.fr;

Contributing authors: drulliere@emse.fr; bay@emse.fr;

Abstract

In the context of Bayesian shape-restricted function estimation, finite-dimensional Gaussian process (GP) approximations, such as those proposed in *Maatouk and Bay* (2017), have gained significant attention due to their flexibility and their ability to handle a wide range of shape constraints. However, these approaches have limitations in higher dimensions, as sampling relies on a regular discretization of the input space. Additionally, their computational complexity grows cubically with the number of observations. In this paper, we propose a new adaptive low-rank methodology that enables finite-dimensional GP approximations to handle multidimensional input spaces, accommodate large datasets, and address multiple and complex shape constraints, even in extreme cases. This approach is based on iteratively adding *optimal*, non-uniform, one-dimensional discretization points and ordering the active input variables. An efficient Markov chain Monte Carlo (MCMC) method for sampling the full posterior distribution, with particular attention paid to large datasets n ≫ 10,000 and hyperparameter selection under shape constraints, is also developed. The proposed method significantly reduces computational complexity and mitigates the *mass-shifting* phenomenon observed in the posterior distribution, thereby substantially improving both prediction accuracy and uncertainty quantification. Its strong empirical performance is demonstrated on several synthetic and real-world datasets.

**Keywords:** Gaussian processes, Complex shape constraints, Multidimensional input, Large datasets, Elliptical slice sampling, MAP estimator, MCMC

# 1 Introduction

Bayesian methods are widely used in signal processing and other research areas due to their robustness in handling noisy data and their efficiency in solving complex problems (Green et al, 2025). In this framework, Gaussian processes (GPs) are known as efficient prior distributions over
function spaces. They are commonly used to address tasks such as regression and classification for
machine learning (Gramacy, 2020; Williams and Rasmussen, 2006). Over time, GPs have become
a popular Bayesian tool across a broad spectrum of applications, with their use expanding into
disciplines such as geosciences, physics, biology, chemistry, engineering, finance, and machine
learning (Blanchet-Scalliet et al, 2023; Chevalier et al, 2014; Finocchio and Schmidt-Hieber,
2023; Kim et al, 2023; Mongwe et al, 2025; Montesinos-López et al, 2023; Petit, 2022). However,
the application of GP models become intractable for large datasets because the time complexity
scales as $O(n^3)$ and the storage requirement as $O(n^2)$, where *n* denotes the number of training
samples. Several efficient methods for handling large datasets have been developed in the
literature, such as variational approaches to sparse GPs (Hensman et al, 2013; Titsias, 2009)
and lower-dimensional subspace methods based on random projections of the dataset (Banerjee
et al, 2013). Moreover, GPs face an additional challenge when incorporating an infinite number
of inequality (shape) constraints (Abrahamsen and Benth, 2001; Kleijnen and Van Beers, 2013),
since the resulting posterior distribution is generally no longer Gaussian. Our objective is to
develop an approach that can simultaneously address both issues, namely the large dataset and
the shape constraints.

Shape constraints such as monotonicity, boundedness and convexity frequently appear in
real-world data applications (Cousin et al, 2016; Johnson and Jiang, 2018; Kochan and Yang,
2025; Zhou et al, 2019). In the context of Bayesian shape-restricted function estimation and
more generally in statistics, incorporating shape constraints significantly enhances prediction
accuracy and yields more realistic confidence intervals (Agrell, 2019; Chipman et al, 2022;
Jidling et al, 2017; Lenk and Choi, 2017; Spiller et al, 2023; Ustyuzhaninov et al, 2020; Riihimaki
and Vehtari, 2010). Advancements in constrained GPs remain an active area of research, with
ongoing efforts to improve their applicability across diverse settings. For instance, the authors
in Lenk et al (2024) have recently developed a Hierarchical Bayesian (HB) model for multi-
group data. Their approach allows for shape constraints and is especially efficient when data
within groups are sparse. Benavoli and Azzimonti (2024) proved that linearly constrained GP
(LCGP) for regression, based on forcing constraints on a finite set of points, can be seen as
Skew Gaussian Processes (SkewGPs).

In the present paper, we consider the finite-dimensional GP approximation developed by
Maatouk and Bay (2017), which has attracted significant attention due to its flexibility and its
ability to incorporate a wide range of linear inequality constraints (shape constraints), such
as monotonicity, convexity, and boundedness. These constraints can be applied individually,
sequentially, or in combination. Unlike some other methods, this approach ensures that the linear
inequality constraints are respected throughout the entire domain, as they are reformulated
as equivalent linear inequality constraints on the basis coefficients. Its theoretical convergence
has been explored in several studies (Bay et al, 2016, 2017; Grammont et al, 2024), within the
context of Reproducing Kernel Hilbert Space (RKHS) (Aronszajn, 1950; Berlinet and Thomas-
Agnan, 2011). Additionally, its performance has been evaluated across various domains, such as
finance (Chataigner et al, 2021; Cousin et al, 2016), nuclear physics (Zhou et al, 2019), nuclear
safety (Maatouk et al, 2025a), and chemistry (Williams et al, 2023). The main idea is based
on a piecewise linear interpolating approximation on a finite-dimensional space using basis
expansions. These basis functions possess attractive properties that are not necessarily shared by
other basis functions such as Bernstein polynomials (Curtis and Ghosh, 2011), regression splines
(Cai and Dunson, 2007; Meyer et al, 2011), and restricted splines (Shively et al, 2011). Through
this approach, the linear inequality constraints (shape constraints) are equivalently reformulated
as finite linear constraints into the basis coefficients. Thus, the resulting posterior distribution
is generated using a truncated multivariate normal (tMVN) distribution (Maatouk and Bay,
2016). However, in high-dimensional input spaces, uniform discretization of the input leads to a
high-dimensional tMVN distribution. Consequently, the finite-dimensional GP approximation
becomes computationally prohibitive and impractical. In this context, several recent works have
been developed (Bachoc et al, 2022; Maatouk et al, 2025a). Nevertheless, these approaches still
face limitations, such as handling complex shape constraints and managing large numbers of
observations in multidimensional settings.

To address this issue, we develop a new methodology that combines adaptive rank reduction
with the ordering of active input variables. More precisely, we first order the active input
variables and iteratively compute *optimal*, non-uniform one-dimensional knots. This significantly
reduces the computational complexity of the GP approximation in multidimensional input
spaces. The proposed methodology is well adapted to the smooth circulant embedding technique
introduced by Ray et al (2020) and further developed by Maatouk et al (2025b, 2024), enabling
it to efficiently handle large datasets n ≫ 10,000 and complex shape constraints. Additionally,
it is particularly well suited for hyperparameter selection processes that incorporate shape
constraints. The contributions of the present paper, compared to previous work, are listed below:

* We investigate the limitations of the finite-dimensional GP approximation in both low- and high-dimensional input spaces when a regular discretization is employed.
* We propose a new methodology based on iteratively computing optimal non-uniform one-dimensional discretization points and ordering active input variables using criteria related to the given data and the maximum a posteriori (MAP) estimator. We additionally develop an efficient MCMC sampler well suited for the proposed approach, which can handle large datasets and complex shape constraints without computing the full posterior conditional distribution.

* We explore the effectiveness of this methodology in several scenarios: multidimensional input settings, large datasets, complex shape constraints in extreme cases, multiple shape constraints, and monotonicity with respect to a subset of the input space.

In Example 1, we empirically demonstrate the necessity of non-uniform discretization for the finite-dimensional GP approximation. In high-dimensional input settings, this approach becomes prohibitive and infeasible when fine uniform discretization is employed. For simplicity of interpretation and visualization, we first consider the one-dimensional case. However, the methodology focuses on multidimensional input settings, and its performance has been demonstrated through several numerical examples.

**Example 1** (Motivating example: one-dimensional input space). Without loss of generality, we focus on monotonicity (non-decreasing) constraints. We consider the underlying function $f: x \mapsto 3/(1 + \exp(-30x + 2.1))$ defined on the unit interval $X$, and observe it under additive Gaussian noise. This is a challenging scenario because the function $f$ grows rapidly and becomes flat on $[0.2, 1]$. Such situations commonly arise when estimating chemical reactions as functions of temperature or when estimating cumulative distribution functions.

In Figure 1, we investigate the impact of the discretization of the input space on the finite-dimensional GP approximation. In the top panels (left and right), a uniform discretization of the input set $X$ is employed with $N = 4$ and $N = 50$ discretization points (knots), respectively (light gray vertical dashed lines). In the bottom panel, an optimal non-uniform discretization is employed using the methodology developed in the present paper. The black stars represent the $n = 100$ training samples generated from (1) using the monotone (non-decreasing) function $f$ (black solid curve) and a zero-mean Gaussian noise $N(0, \sigma^2)$, where $\sigma = 0.5$. The covariates $\{x_i\}$, $i = 1, \dots, n$ are generated using the Latin-hypercube sampling (LHS) design from `lhsDesign` function implemented in the R-package `DiceDesign` (Dupuy et al., 2015). The two red and blue (dash-dotted and dash) curves represent the mAP and MAP estimator, respectively (Definitions 1 and 2 below). The light gray vertical dashed lines represent the locations of the input discretization points. The gray shaded area represents the 95% confidence interval obtained by
<figure><img src="image_3.png" /><figcaption>Fig. 1 GP approximation under monotonicity constraints with uniform discretization (top left and right, with N = 4 and N = 50, respectively) and optimal non-uniform discretization (bottom).</figcaption></figure>

generating monotonicity sample paths using 6,000 MCMC iterations, with the first 1,000 iterations discarded as burn-in. The MSE\_mAP and MSE\_MAP represent, respectively, the mean squared error (MSE) criterion when using the mAP and MAP estimators:

$$MSE = \frac{1}{n_{test}} \sum_{i=1}^{n_{test}} (f(x_{i}^{test}) - \hat{y}_{i})^{2},$$

where $n_{test} = 1000$, $\{x_{i}^{test}\}$ are equally spaced test points, and $\hat{y}_{i}$ is the predictive value at $x_{i}^{test}$. First, with only a few uniform knots (top-left panel), the finite-dimensional GP approximation provides a poor predictive model. Second, we observe that fine uniform discretization (top-right panel) leads to an overfitting and accentuates the mass-shifting phenomenon of the posterior distribution and leads to a biased Bayesian estimator, especially in the flat regions of the underlying function. It is important to mention that the mass-shifting phenomenon occurs in the low-dimensional marginal densities of a dependent zero-mean multivariate normal (MVN) distribution restricted to the positive orthant, as well quantified by Zhou et al (2024). It is worth noting that the MAP (red dashed curve) is not impacted by this phenomenon, as this estimator relies only on a quadratic optimization problem (Maatouk et al, 2025b, Sect. 3.2). Third, we observe that non-uniform discretization significantly improves both prediction accuracy and coverage intervals, and a small number of optimal discretization points are sufficient to achieve a good approximation of the underlying function.
The article is structured as follows: in Section 2 we provide a brief review of GP regression (GPR). In Section 3, we develop a finite-dimensional GP approximation under shape constraints and noisy observations for multidimensional input settings, using a non-uniform discretization. We also define the associated posterior estimators and analyze their corresponding computational complexities. In Section 4, we present a grid-based strategy for low-rank representations and for ordering active input variables. Furthermore, we develop an efficient MCMC sampler based on the smooth circulant embedding technique for generating the full posterior distribution. The performance of the proposed approach is demonstrated using both synthetic data in Section 5 and real-world data studies in Section 6.

## 2 Gaussian process regression (GPR) review

In this section, we provide a brief overview of the GPR framework, a probabilistic approach to nonparametric function estimation. Specifically, we consider the following regression model:

$$y_i = f(\mathbf{x}_i) + \epsilon_i. \qquad (1)$$

where $f: \mathcal{X} \subset \mathbb{R}^d \to \mathbb{R}$ denotes an unknown real-valued function that generates the observations $\mathbf{y} = [y_1, \dots, y_n]^\top \in \mathbb{R}^n$, and $\mathbf{x}_i \in \mathcal{X}$ is a covariate of dimension $d \ge 1$. In practice, the observed function values are often noisy, and it is typically assumed that the observation noise is additive (1), where the noise terms $\epsilon_i \stackrel{\text{i.i.d.}}{\sim} \mathcal{N}(0, \sigma^2)$, $i = 1, \dots, n$ are assumed to be independent and identically distributed (i.i.d.) Gaussian random variables with zero-mean and constant variance of $\sigma^2$. The reader is referred to Binois et al (2018); Binois and Gramacy (2021) for a detailed discussion of input-dependent noise variance. The set $\mathcal{X}$ is a compact subset of $\mathbb{R}^d$ and without loss of generality, we suppose that $\mathcal{X}$ is a unit hypercube. We denote by $\{(\mathbf{x}_i, y_i)\}_{i=1}^n$ the set of $n$ training samples and by $\mathbb{X}$ the $n \times d$ design matrix obtained by concatenating the covariates, i.e., $\mathbb{X} = [\mathbf{x}_1, \dots, \mathbf{x}_n]^\top$. In matrix form, the regression problem (1) can be expressed as:

$$\mathbf{y} = f(\mathbb{X}) + \mathbf{\epsilon},$$

where $\mathbf{\epsilon} = [\epsilon_1, \dots, \epsilon_n]^\top$ is a zero-mean Gaussian noise vector with covariance matrix $\sigma^2 I_n$, with $I_n$ the $n \times n$ identity matrix.

Let $(Z(\mathbf{x}))_{\mathbf{x} \in \mathcal{X}}$ be a stochastic process on $\mathcal{X}$ such that $E[Z(\mathbf{x})^2] < +\infty$, for all $\mathbf{x} \in \mathcal{X}$. Then, the function $\mu: \mathcal{X} \to \mathbb{R}$, with $\mu(\mathbf{x}) = E[Z(\mathbf{x})]$ is called the mean function of $Z$ and the
function $k: \mathcal{X} \times \mathcal{X} \rightarrow \mathbb{R}$, with $k(\mathbf{x}, \mathbf{x}') = \text{Cov}(Z(\mathbf{x}), Z(\mathbf{x}'))$ is called the covariance function of $Z$.

GPR is a nonparametric and probabilistic approach to regression. It is based on placing a GP prior over the unknown function $f$ in (1) (Williams and Rasmussen, 2006, Section 2.2). A GP is a stochastic process (a collection of random variables) indexed by time or space, where every finite subset of these variables follows a MVN distribution. Suppose that $(Z(\mathbf{x}))_{x \in \mathcal{X}}$ is a GP with mean function $\mu(\cdot)$ and covariance function $k(\cdot, \cdot)$, i.e., $Z \sim \mathcal{GP}(\mu(\cdot), k(\cdot, \cdot))$. Thus,

$$Z(\mathbf{x}) = \mu(\mathbf{x}) + Y(\mathbf{x}), \quad \forall \mathbf{x} \in \mathcal{X},$$

where $Y$ is a zero-mean GP with covariance function $k(\cdot, \cdot)$, i.e., $Y \sim \mathcal{GP}(0, k(\cdot, \cdot))$. In the GPR framework, the selection of the positive-definite and bounded covariance function (or kernel) $k$ is critical. This function measures the similarity between two data points. Additionally, it encodes assumptions about the underlying process, such as differentiability and periodicity, through its properties. The covariance kernel is typically selected from a stationary parametric family such as Matérn family (Genton, 2001). A stationary kernel can be expressed as a function of the difference between two input points, $\mathbf{x} - \mathbf{x}'$, whereas a nonstationary kernel depends on the individual values of $\mathbf{x}$ and $\mathbf{x}'$ themselves, rather than solely on their difference. In this paper, we consider a tensor product of one-dimensional Matérn covariance kernels with smoothness parameter $\nu$, defined by

$$k(x, x') = \tau^2 \frac{2^{1-\nu}}{\Gamma(\nu)} \left( \frac{\sqrt{2\nu}}{\ell} |x - x'| \right)^{\nu} B_{\nu} \left( \frac{\sqrt{2\nu}}{\ell} |x - x'| \right), \quad (2)$$

for any $x, x' \in \mathcal{X} = [0, 1]$, where $\Gamma(\cdot)$ is the Gamma function and $B_{\nu}(\cdot)$ denotes the modified Bessel function of the second kind of order $\nu$. The positive parameter $\nu$ regulates the degree of smoothness of the GP sample paths. A process with the Matérn kernel of order $\nu$ is $[\nu - 1]$ times differentiable (Williams and Rasmussen, 2006, Sect. 4.2.1). The positive kernel parameters $\tau^2$ and $\ell$ are referred to as the *signal variance* and the *correlation length-scale*, respectively. They are often estimated using maximum likelihood or cross-validation techniques (see, e.g., Bachoc (2013) and (Williams and Rasmussen, 2006, Sect. 5.4)). In this paper, we develop an efficient strategy for hyperparameter selection that accounts for shape constraints.
## 2.1 Unconstrained Gaussian process regression

In this section, we consider only noisy observations, excluding the shape constraints. Conditioned on the observations $\mathbf{y} = [y_1, \dots, y_n]^\top$, the resulting process remains a GP (Williams and Rasmussen, 2006, Sect. 2.2), i.e.,

$$ \{Y(\cdot)|\mathbf{y}\} \sim \mathcal{GP}(\tilde{\mu}(\cdot), \tilde{k}(\cdot, \cdot)), $$

where the conditional mean function $\tilde{\mu}$ and covariance function $\tilde{k}$ at any new location $\mathbf{x} \in \mathcal{X}$ are given by:

$$ \left\{ \begin{array}{l} \tilde{\mu}(\mathbf{x}) = \mathrm{E}[Y(\mathbf{x})|\mathbf{y}] = k(\mathbf{x}, \mathbf{X})^{\top} (k(\mathbf{X}, \mathbf{X}) + \sigma^2 \mathrm{I}_n)^{-1} \mathbf{y}; \\ \tilde{k}(\mathbf{x}, \mathbf{x}') = k(\mathbf{x}, \mathbf{x}') - k(\mathbf{x}, \mathbf{X})^{\top} (k(\mathbf{X}, \mathbf{X}) + \sigma^2 \mathrm{I}_n)^{-1} k(\mathbf{x}', \mathbf{X}); \end{array} \right. \qquad (3) $$

with $\mathrm{I}_n$ the $n \times n$ identity matrix. Let us recall that $k(\mathbf{X}, \mathbf{X})$ is the $n \times n$ covariance matrix of the Gaussian vector $Y(\mathbf{X}) \in \mathbb{R}^n$, and $k(\mathbf{x}, \mathbf{X}) = [k(\mathbf{x}, \mathbf{x}_1), \dots, k(\mathbf{x}, \mathbf{x}_n)]^\top$ is the vector of covariance between $Y(\mathbf{x}) \in \mathbb{R}$ and $Y(\mathbf{X}) \in \mathbb{R}^n$. Here, we denote by $k(\mathcal{E}, \mathcal{F})$ the matrix whose entries correspond to the covariance function evaluated at all pairs of points, with one point from the set $\mathcal{E}$ and the other from the set $\mathcal{F}$. The posterior mean function $\tilde{\mu}(\cdot)$ in (3) is commonly employed as the surrogate model, while the posterior covariance function $\tilde{k}(\cdot, \cdot)$ quantifies the uncertainty associated with the surrogate's predictions. It is worth noting that noise-free GPR, where $f(\mathbf{x}_i) = y_i$ can be considered as an interpolation method, as it exactly fits all observed data points. In contrast, GPR with observation noise, as described in (1), permits deviations from the data, thereby accommodating fitting errors. Finally, the equations in (3) are commonly referred to as the *kriging equations* in geostatistics. The posterior mean $\tilde{\mu}$ describes the best (minimizing mean-squared prediction error MSPE) linear unbiased predictor (BLUP).

## 2.2 Constrained Gaussian process regression

In the following, we consider Bayesian shape-restricted function estimation using a GP prior. This Bayesian approach is known as constrained GPR or linearly restricted GPR (Benavoli and Azzimonti, 2024; Jidling et al, 2017) which is widely investigated in the machine learning community. The underlying function is assumed to satisfy shape constraints, such as monotonicity, convexity, and boundedness, as well as any combination of these (referred to as multiple shape constraints). Throughout this paper, we denote by $\mathcal{C}$ the set of functions satisfying shape
constraints. Thus, we are interested by studying the following posterior distribution

$$ \{Y(\cdot) | \mathbf{y}, Y \in \mathcal{C}\}. \qquad (4) $$

A common approach to incorporate shape constraints into a GP model is to employ a Bayesian linear model with a suitable basis expansion, such as those mentioned in the introduction. In this work, we adopt the finite-dimensional GP approximation proposed by Maatouk and Bay (2017), which employs a different type of basis function. The extension to multidimensional input settings is achieved via a tensor product of hat basis functions. We slightly modify these basis functions from a symmetric to an asymmetric form to adapt the approach to a non-uniform subdivision. Importantly, the methodology developed here is flexible and can be extended to accommodate other types of shape constraints and basis functions, such as those introduced by Maatouk et al (2025a).

# 3 Gaussian process approximations

In this section, we consider the finite-dimensional GP approximation proposed by Maatouk and Bay (2017), which is capable of incorporating a wide range of shape constraints over the entire domain. This approach can be interpreted as a Bayesian linear model, often referred to as the weight-space view of GPs (Williams and Rasmussen, 2006, Sect. 2.1.1). Subsection 3.1 reviews the one-dimensional case with uniform discretization, while Subsection 3.2 extends existing methods to the multidimensional case with non-uniform discretization.

## 3.1 One-dimensional case with uniform discretization

In this section, we consider the one-dimensional case, i.e., $d = 1$, where $\mathcal{X}$ is the unit interval. We let $s$ denote a sequence of $N \ge 2$ equally spaced knots $t_j \in \mathcal{X}$, for $j = 1, \dots, N$, with spacing $\delta_N = 1/(N-1)$. Thus, $s = (t_1, \dots, t_N)$, where $t_1 = 0$ and $t_N = 1$. Consequently, $t_j = (j-1)\delta_N$ and $t_{j+1} - t_j = \delta_N$. To simplify the notation while preserving an equally spaced subdivision, we set $t_0 = -\delta_N$ and $t_{N+1} = 1 + \delta_N$. Additionally, we denote by $\phi_j$ the one-dimensional symmetric hat basis function corresponding to the knot $t_j$, which has the following form:

$$ \phi_j(x) = \begin{cases} \frac{x-t_{j-1}}{\delta_N} & \text{if } x \in [t_{j-1}, t_j]; \\ \frac{t_{j+1}-x}{\delta_N} & \text{if } x \in [t_j, t_{j+1}]; \\ 0 & \text{otherwise;} \end{cases} $$
for any $j = 1, \dots, N$, i.e., the B-spline of degree 1 corresponding to the knots $\{t_{j-1}, t_j, t_{j+1}\}$ (Eilers and Marx, 1996, Fig. 1.a). It is worth noting that for any $i, j = 1, \dots, N$, we have $\phi_j(t_i) = \delta_{i,j}$, where $\delta_{i,j}$ is the Kronecker delta, equal to one if $i = j$ and zero otherwise. Furthermore, $\sum_{j=1}^N \phi_j(x) = 1$, for any $x \in \mathcal{X}$. These two properties are used to demonstrate the convergence of the proposed GP approximation.

The finite-dimensional GP approximation is inspired by the following piecewise linear interpolation theory.

**Lemma 1** (Uniform convergence of piecewise linear interpolation). Let $f$ be a continuous function on $\mathcal{X}$, i.e., $f \in C^0(\mathcal{X}, \mathbb{R})$. Then, the piecewise linear interpolating function $\tilde{f}_N(\cdot) = \sum_{j=1}^N f(t_j)\phi_j(\cdot)$ converges uniformly to $f$, when $N$ tends to infinity. Here, the sequence $\{t_j\}$ represents the $N$ equally spaced knots of the input set $\mathcal{X}$, where $t_j = (j-1)/(N-1)$, for any $j = 1, \dots, N$. Furthermore, the functions $\tilde{f}_N$ and $f$ coincide at each knot $t_j$, i.e., $\tilde{f}_N(t_j) = f(t_j)$, which implies that the approximation satisfies the interpolation property (see the black points in Figure 2 for $N = 5$).

<figure><img src="image_4.png" /><figcaption>**Fig. 2** Approximating a function $f$ (red solid curve) using a piecewise linear interpolating function $\tilde{f}_N$ (black dashed curve). A uniform subdivision is used with $N = 5$ hat basis functions (gray dash-dotted triangles).</figcaption></figure>

Based on Lemma 1, expanding the finite-dimensional GP approximation $\tilde{Y}_N$ of $Y$ in the interpolation basis implies the model

$$ \tilde{Y}_N(x) = \sum_{j=1}^{N} Y(t_j)\phi_j(x) = \phi(x)^{\top}\xi, \quad x \in \mathcal{X}, \qquad (5) $$

where $\xi = [\xi_1, \dots, \xi_N]^\top$ with $\xi_j = Y(t_j)$, and $\phi(\cdot) = [\phi_1(\cdot), \dots, \phi_N(\cdot)]^\top \in \mathbb{R}^N$. Since the parent process $Y$ is assumed to be a zero-mean Gaussian with covariance function $k(\cdot, \cdot)$, then $\xi$ is a zero-mean Gaussian vector with covariance matrix $\tau^2 K$, i.e., $\xi \sim N(0, \tau^2 K)$, where
0 = [0, ..., 0]ᵀ is the N-dimensional zero vector and

$$ \tau^2 \mathbf{K}_{i,j} = \mathrm{Cov}(\xi_i, \xi_j) = \mathrm{Cov}(Y(t_i), Y(t_j)) = k(t_i, t_j), \quad i, j = 1, \dots, N, \qquad (6) $$

with $\tau^2 > 0$ denoting the signal variance, which controls the variance of the prior. Moreover, for any $x, x' \in \mathcal{X}$, we have

$$ k^N(x, x') = \mathrm{Cov}(\tilde{Y}_N(x), \tilde{Y}_N(x')) = \tau^2 \phi(x)^\top \mathbf{K} \phi(x'), \qquad (7) $$

which represents the covariance function of $\tilde{Y}_N$. The matrix $\mathbf{K}$ has a rank $N$, usually smaller than the dimension of the dataset, which motivates referring to $\tilde{Y}_N$ as a low-rank representation of the process $Y$.

**Lemma 2.** If the covariance function $k$ of the parent GP $Y$ is positive-definite, then the covariance function $k^N$ of the finite-dimensional GP defined in (7) is also positive-definite.

Unlike some other Bayesian linear approaches, such as those mentioned in the introduction, the Gaussian weights or basis coefficients {$\xi_j$} in (5) are correlated, resulting in greater flexibly. Moreover, their correlation is derived from the correlation function of the parent GP $Y$ (6). $\tilde{Y}_N$ is a finite-dimensional GP, and converges uniformly to $Y$ as $N$ tends to infinity, with probability one. This implies that every realization of $\tilde{Y}_N$ converges uniformly to a realization of $Y$. From (6), if a stationary covariance function $k(\cdot, \cdot)$, such as one from the Matérn family (Genton, 2001), is employed (i.e., $k(t, t') = k(|t - t'|)$, $\forall t, t' \in \mathbb{R}$), then equally spaced knots {$t_j$} in the input domain $\mathcal{X}$ result in a Toeplitz covariance matrix $\mathbf{K}$. This idea has been explored in both Ray et al (2020) and Maatouk et al (2025c), where the prior was generated efficiently using the Fast Fourier Transform (FFT) (Wood and Chan, 1994) and fast large-scale approaches developed by Maatouk et al (2023), respectively.

**Lemma 3.** Assume that $Y$ is a GP with mean function $\mu(\cdot)$ and covariance function $k(\cdot, \cdot)$. Suppose also that $\mu(\cdot)$ is at least a continuous function on $\mathcal{X}$. Then, $\tilde{\mu}_N(\cdot) := E[\tilde{Y}_N(\cdot)]$ converges uniformly to $\mu(\cdot)$, where $\tilde{Y}_N$ is the finite-dimensional GP approximation defined in (5). Additionally, if $\mu$ is Lipchitz, then $\|\tilde{\mu}_N - \mu\|_{\infty} \le M/(N - 1)$, where $M$ is a finite positive constant such that $|\mu(x) - \mu(y)| \le M|x - y|$, for any $x, y \in \mathcal{X}$.

As shown in Figure 1 (top panels), using a uniform discretization leads to a biased model (top-right panel) when the number of knots $N$ is too high due to the mass-shifting phenomenon and to a poor predictive model (top-left panel) when only a few uniform knots are used, i.e.,
when N is low.

To address these issues, we develop in the following sections the finite-dimensional GP approximation (5) in multidimensional case with optimal non-uniform discretization.

## 3.2 Multidimensional case with non-uniform discretization

In this section, we extend existing uniform grid of Maatouk and Bay (2017) to non-uniform discretizations in the multidimensional setting. The input set $\mathcal{X}$ is the unit hypercube, i.e., $\mathcal{X} = [0,1]^d$ with $d \ge 2$. We consider the finite-dimensional GP approximation (5) in the multidimensional case with a non-uniform discretization of $\mathcal{X}$. This leads to the following model

$$ 
\begin{aligned}
\tilde{Y}_S(\mathbf{x}) &:= \sum_{j_1=1}^{N_1} \dots \sum_{j_d=1}^{N_d} Y(t_{j_1}^1, \dots, t_{j_d}^d) [\phi_{j_1}^1(x_1) \times \dots \times \phi_{j_d}^d(x_d)], \\
&= \sum_{j_1=1}^{N_1} \dots \sum_{j_d=1}^{N_d} \xi_{j_1, \dots, j_d} [\phi_{j_1}^1(x_1) \times \dots \times \phi_{j_d}^d(x_d)], \\
&= \Phi_S(\mathbf{x})^\top \xi_S, \quad \mathbf{x} = (x_1, \dots, x_d) \in \mathcal{X},
\end{aligned}
\quad (8) $$

where $\xi_S = \{\xi_{j_1, \dots, j_d}\}$ is a $\prod_{\zeta=1}^d N_\zeta$-dimensional vector, with $d \ge 2$, and $S = (s^{(1)}, \dots, s^{(d)})$ is the subdivision of the input set $\mathcal{X}$. For example, for any $\zeta \in \{1, \dots, d\}$, we have $s^{(\zeta)} = \{t_1^\zeta, \dots, t_{N_\zeta}^\zeta\}$, representing a given ordered non-uniform partition of the input dimension $\zeta$, i.e., $t_1^\zeta < \dots < t_{N_\zeta}^\zeta$, with $N_\zeta \ge 2$. Clearly, the proposed GP approximation $\tilde{Y}_S$ defined in (8) depends on the subdivision $S$, as the parent GP $Y$ is evaluated on $S$ and the basis functions $\phi_{j_\zeta}^\zeta$ themselves depend on it. Accordingly, we include $S$ as subscript in $\tilde{Y}_S$, $\Phi_S$ and $\xi_S$ to explicitly indicate their dependence on the subdivision. The originality of the model (8) lies in the fact that the Gaussian weight vector $\xi_S$ is derived from the parent GP $Y$ via the relation $\xi_S := Y(t_{j_1}^1, \dots, t_{j_d}^d)$, where $(t_{j_1}^1, \dots, t_{j_d}^d) \in \mathcal{X}$ denotes a $d$-dimensional knot corresponding to a given multi-index $j_1, \dots, j_d$. Therefore, the $\prod_{\zeta=1}^d N_\zeta$ Gaussian weights $\{\xi_{j_1, \dots, j_d}\}$ are correlated. Here, $\phi_{j_\zeta}^\zeta(\cdot)$ denotes a one-dimensional asymmetric hat basis function (9) corresponding to the input dimension $\zeta$ and the one-dimensional knots $\{t_{j_\zeta-1}^\zeta, t_{j_\zeta}^\zeta, t_{j_\zeta+1}^\zeta\}$, where $\zeta \in \{1, \dots, d\}$. Additionally, $N_\zeta \ge 2$ represents the number of one-dimensional knots for the input dimension $\zeta$ or the size of the one-dimensional subdivision $s^\zeta = \{t_1^\zeta, \dots, t_{N_\zeta}^\zeta\}$. We later explain how optimal non-uniform knots, $\{t_{j_\zeta}^\zeta\}$, are computed, for any $j_\zeta \in \{1, \dots, N_\zeta\}$ and $\zeta \in \{1, \dots, d\}$. Throughout this paper, superscripts and subscripts in the variables $t$ and $\phi$ indicate the input dimension and the corresponding single-index, respectively. Let $t_0^\zeta < 0$ and $t_{N_\zeta+1}^\zeta > 1$ be two arbitrary values
outside the unit interval $[0, 1]$, which have no effect on the model itself. Their role is solely to adjust the definition of the asymmetric one-dimensional *hat* basis functions at the boundaries, as described below. The one-dimensional asymmetric *hat* basis function $\phi_{j\zeta}^{\zeta}$ associated to the one-dimensional knots $\{t_{j\zeta-1}^{\zeta}, t_{j\zeta}^{\zeta}, t_{j\zeta+1}^{\zeta}\}$ is given by

$$ \phi_{j\zeta}^{\zeta}(x) = \begin{cases} \frac{x-t_{j\zeta-1}^{\zeta}}{t_{j\zeta}^{\zeta}-t_{j\zeta-1}^{\zeta}} & \text{if } x \in [t_{j\zeta-1}^{\zeta}, t_{j\zeta}^{\zeta}]; \\ \frac{t_{j\zeta+1}^{\zeta}-x}{t_{j\zeta+1}^{\zeta}-t_{j\zeta}^{\zeta}} & \text{if } x \in [t_{j\zeta}^{\zeta}, t_{j\zeta+1}^{\zeta}]; \\ 0 & \text{otherwise}; \end{cases} \quad (9) $$

for any $j_\zeta = 1, \dots, N_\zeta$ and any $\zeta \in \{1, \dots, d\}$, with $N_\zeta \ge 2$. Clearly, a function of the form $\phi_{j\zeta}^{\zeta}$ is a *hat* basis function centered at $t_{j\zeta}^{\zeta}$, with compact support $[t_{j\zeta-1}^{\zeta}, t_{j\zeta+1}^{\zeta}]$. It is asymmetric if $t_{j\zeta+1}^{\zeta} - t_{j\zeta}^{\zeta} \ne t_{j\zeta}^{\zeta} - t_{j\zeta-1}^{\zeta}$. Figure 3 shows the case when $N_\zeta = 4$ and $s^\zeta = \{0, 0.25, 0.75, 1\}$. The gray dash-dotted triangles represent the basis functions $\{\phi_{j\zeta}^{\zeta}\}$, for the single-index $j_\zeta = 1, \dots, 4$. For any $\zeta \in \{1, \dots, d\}$, we have

$$ t_{j\zeta}^{\zeta} = \begin{cases} 0 & \text{if } j_\zeta = 1; \\ 1 & \text{if } j_\zeta = N_\zeta. \end{cases} $$

Furthermore, $\phi_{j\zeta=j}^{\zeta}(t_{j\zeta}^{\zeta}=i) = \delta_{i,j}$ for any $i, j \in \{1, \dots, N_\zeta\}$, where $\delta_{i,j}$ is the Kronecker delta, equal to one if $i = j$ and zero otherwise. It is worth noting that the non-uniform discretization of the input set $\mathcal{X}$ preserves the following properties:

**Lemma 4.** For any $N_\zeta \ge 2$, where $\zeta \in \{1, \dots, d\}$, we have

$$ \sum_{j_\zeta=1}^{N_\zeta} \phi_{j_\zeta}^{\zeta}(x_\zeta) = 1, \quad \forall x_\zeta \in [t_1^\zeta, t_{N_\zeta}^\zeta], \quad (10) $$

which is known as partition of unity property. Moreover, we have

$$ \sum_{j_1=1}^{N_1} \cdots \sum_{j_d=1}^{N_d} \left( \prod_{\zeta=1}^{d} \phi_{j_\zeta}^{\zeta}(x_\zeta) \right) = 1, \quad (11) $$

for all $x = (x_1, \dots, x_d) \in \mathcal{X}$, with $d \ge 2$. Clearly, for a given multi-index $j_1, \dots, j_d$, the tensor product $\prod_{\zeta=1}^d \phi_{j_\zeta}^{\zeta}(\cdot)$ is equal to one at the $d$-dimensional knot $(t_{j_1}^1, \dots, t_{j_d}^d)$, and its compact support is the hypercube $\prod_{\zeta=1}^d [t_{j_\zeta-1}^{\zeta}, t_{j_\zeta+1}^{\zeta}]$.

Figure 3 illustrates the approximation of a given continuous function $f$ (red solid curve) using $\tilde{f}_{N_1}(\cdot) = \sum_{j_1=1}^{N_1} f(t_{j_1}^1) \phi_{j_1}^1(\cdot)$ (black dashed curve), as described in Lemma 1. This is analogous to the approximation shown in Figure 2, but here the discretization is non-uniform, with $N_1 = 4$. We removed only the middle knot at $t_3 = 0.5$, since the target function is approximately linear
<figure><img src="image_5.png" /><figcaption>Fig. 3 Same settings as in Figure 2, but using a non-uniform subdivision with $N_1 = 4$ hat basis functions (gray dash-dotted triangles).</figcaption></figure>

between the second knot $t_2 = 0.2$ and the fourth knot $t_4 = 0.75$. Thus, the new non-uniform subdivision is $t_1^1 = 0$, $t_2^1 = 0.25$, $t_3^1 = 0.75$, and $t_4^1 = 1$. Notably, the prediction accuracy remains nearly unaffected by the removal of a single knot.

**Lemma 5.** Assume that $f$ is a Lipschitz function, i.e., there exists a finite constant $M > 0$ such that $|f(x) - f(y)| \le M|x - y|$, for any $x, y \in \mathcal{X}$. Let $\tilde{f}_{N_1}(x) = \sum_{j_1=1}^{N_1} f(t_{j_1}^1)\phi_{j_1}^1(x)$, $x \in \mathcal{X}$, where $N_1 \ge 2$, $\{t_{j_1}^1\}$ is a sequence of non-uniform one-dimensional knots, and $\phi_{j_1}^1(\cdot)$ is the one-dimensional basis function corresponding to the non-uniform knots $\{t_{j_1-1}^1, t_{j_1}^1, t_{j_1+1}^1\}$. Suppose also that $\max_{j_1}(t_{j_1+1}^1 - t_{j_1}^1) \le \kappa/(N_1 - 1)$, $j_1 = 1, \dots, N_1 - 1$, with $\kappa + 1 \le N_1$. Then, $\|f - \tilde{f}_{N_1}\|_{\infty} \le M\kappa/(N_1 - 1)$.

**Corollary 1.** We consider the function approximation

$$ \tilde{f}_{\mathcal{S}}(\mathbf{x}) = \sum_{j_1=1}^{N_1} \cdots \sum_{j_d=1}^{N_d} f(t_{j_1}^1, \dots, t_{j_d}^d) [\phi_{j_1}^1(x_1) \times \cdots \times \phi_{j_d}^d(x_d)], $$

for $\mathbf{x} = (x_1, \dots, x_d) \in \mathcal{X}$, where each $\phi_{j_\zeta}^\zeta$ is an asymmetric hat basis function defined on a grid $\{t_{j_\zeta}^\zeta\}_{j_\zeta=1}^{N_\zeta}$ over the unit interval, and the function $f$ is continuous on the domain $\mathcal{X} = [0, 1]^d$. Suppose that $f$ verifies the following assumptions:

(i) $f$ is uniformly continuous on the compact domain $\mathcal{X}$.

(ii) Each family $\{\phi_{j_\zeta}^\zeta\}_{j_\zeta=1}^{N_\zeta}$ satisfies the partition of unity property:

$$ \sum_{j_\zeta=1}^{N_\zeta} \phi_{j_\zeta}^\zeta(x_\zeta) = 1, \quad \forall x_\zeta \in [0, 1], $$

and each basis function $\phi_{j_\zeta}^\zeta$ has local support.
(iii) The grids {$t_{j_\zeta}^\zeta$} are quasi-uniform: there exists a constant $C > 0$ such that

$$h_{\zeta} := \max_{j_{\zeta}}(t_{j_{\zeta}+1}^{\zeta} - t_{j_{\zeta}}^{\zeta}) \le C \min_{j_{\zeta}}(t_{j_{\zeta}+1}^{\zeta} - t_{j_{\zeta}}^{\zeta}), \quad \forall \zeta = 1, \dots, d.$$

Then, $\tilde{f}_S \to f$ uniformly on $\mathcal{X}$. Additionally, if $f \in H^\alpha(\mathcal{X})$, with $0 < \alpha \le 1$, then the approximation $\tilde{f}_S$ converges weakly to $f$ in $L^2$ as $h \to 0$, with rate $O(h^\alpha)$, where $h = \max_\zeta h_\zeta$, and $H^\alpha(\mathcal{X})$ is the Hölder-Sobolev space.

The proof is deferred to the Appendix. In the next section, we explain how Model (8) can deal with noisy observations and shape constraints. Later in Section 4, we explain how optimal non-uniform knots {$t_{j_\zeta}^\zeta$} are computed, for any $j_\zeta \in \{1, \dots, N_\zeta\}$ and $\zeta \in \{1, \dots, d\}$, and how this model accommodates large datasets, complex shape constraints, and multidimensional input settings.

## 3.3 Constrained GP approximation with noisy observations

In this section, we assume that the GP approximation $\tilde{Y}_S$ defined in (8) satisfies noisy observations, i.e., $\{\tilde{Y}_S(\mathcal{X}) + \epsilon = \mathbf{y}\}$, and shape constraints, i.e., $\tilde{Y}_S(\mathbf{x}) \in \mathcal{C}$. We first focus on the noisy observations. The matrix $\mathcal{X}$ is obtained by concatenating the covariates $\mathbf{x}_i \in \mathbb{R}^d$, i.e., $\mathcal{X} = [\mathbf{x}_1, \dots, \mathbf{x}_n]^\top$. According to (8), they can be expressed in matrix form as follows:

$$\mathbf{X}_S \xi_S + \epsilon = \mathbf{y},$$

where $\mathbf{X}_S$ is the $n \times \prod_{\zeta=1}^d N_\zeta$ design matrix defined by $\mathbf{X}_S := \Phi_S(\mathcal{X})^\top$, and $\epsilon = [\epsilon_1, \dots, \epsilon_n]^\top$ is the noise Gaussian vector. Here, $\prod_{\zeta=1}^d N_\zeta$ represents the total number of non-uniform (one-dimensional) knots, and $n$ denotes the number of observations. We are interested in the situation where $n \gg \prod_{\zeta=1}^d N_\zeta$. The conditional distribution $\{\tilde{Y}_S|\mathbf{X}_S, \mathbf{y}\}$ can be obtained from the conditional distribution on the Gaussian weight, i.e., $\{\xi_S|\mathbf{X}_S, \mathbf{y}\}$:

$$\{\tilde{Y}_S(\cdot)|\mathbf{X}_S, \mathbf{y}\} \stackrel{d}{=} \Phi_S(\cdot)^\top \{\xi_S|\mathbf{X}_S, \mathbf{y}\},$$

where $\xi_S$ is the $\prod_{\zeta=1}^d N_\zeta$-dimensional Gaussian vector with mean zero and covariance matrix $\tau^2 \mathbf{K}_S$, where

$$\tau^2 \mathbf{K}_S = (\mathrm{Cov}(Y(t_{i_1}^1, \dots, t_{i_d}^d), Y(t_{j_1}^1, \dots, t_{j_d}^d))) = k((t_{i_1}^1, \dots, t_{i_d}^d), (t_{j_1}^1, \dots, t_{j_d}^d)),$$
for any $i_ζ$ and $j_ζ$ in $\{1, \dots, N_ζ\}$ with $ζ \in \{1, \dots, d\}$.

On the one hand, we have the following well known result:

**Lemma 6.** Let $\xi_S$ be a $\prod_{\zeta=1}^d N_ζ$-dimensional Gaussian vector with mean zero and covariance matrix $\tau^2 K_S$, where $\tau^2$ is a positive kernel parameter. Suppose also that $\mathbf{y}$ is a given $n$-dimensional vector representing data such that $\mathbf{y} = \mathbf{X}_S\xi_S + \boldsymbol{\epsilon}$, with $\mathbf{X}_S$ an $n \times \prod_{\zeta=1}^d N_ζ$ design matrix and $\boldsymbol{\epsilon} \sim \mathcal{N}(\mathbf{0}, \sigma^2\mathbf{I}_n)$. Then

$$ \{\xi_S | \mathbf{X}_S, \mathbf{y}\} \sim \mathcal{N}(\tilde{\boldsymbol{\mu}}_S, \tilde{\boldsymbol{\Sigma}}_S), \quad \text{where,} $$

$$ \left\{ \begin{array}{l} \tilde{\boldsymbol{\mu}}_S = \tau^2 (\mathbf{X}_S \mathbf{K}_S)^\top (\tau^2 \mathbf{X}_S \mathbf{K}_S \mathbf{X}_S^\top + \sigma^2 \mathbf{I}_n)^{-1} \mathbf{y}; \\ \tilde{\boldsymbol{\Sigma}}_S^{-1} = [\tau^2 \mathbf{K}_S - \tau^4 (\mathbf{X}_S \mathbf{K}_S)^\top (\tau^2 \mathbf{X}_S \mathbf{K}_S \mathbf{X}_S^\top + \sigma^2 \mathbf{I}_n)^{-1} (\mathbf{X}_S \mathbf{K}_S)]^{-1}; \end{array} \right. \quad (12) $$

with $\tilde{\boldsymbol{\mu}}_S$ and $\tilde{\boldsymbol{\Sigma}}_S$ the posterior mean and covariance matrix, respectively.

The conditional distribution in (12) is closely related to the posterior predictive formulas of the conditional GP in (3). It is worth noting that the posterior covariance matrix $\tilde{\boldsymbol{\Sigma}}_S$ is not Toeplitz/block-Toeplitz in general, even when the prior covariance matrix $\tau^2 K_S$ is Toeplitz/block-Toeplitz.

On the other hand, and according to (Williams and Rasmussen, 2006, Sect. 2.1.1) and Bayes' rule, the posterior distribution of the Gaussian weight $\xi_S$ is given by

$$ p(\xi_S | \mathbf{X}_S, \mathbf{y}) = \frac{p(\mathbf{y} | \mathbf{X}_S, \xi_S) p(\xi_S)}{p(\mathbf{y} | \mathbf{X}_S)}, \quad (13) $$

where $p(\mathbf{y} | \mathbf{X}_S)$ is the normalizing constant, also known as the marginal likelihood. It is independent of the weight $\xi_S$ and given by

$$ p(\mathbf{y} | \mathbf{X}_S) = \int_{\mathbb{R}} \prod_{\zeta=1}^d N_\zeta p(\mathbf{y} | \mathbf{X}_S, \xi_S) p(\xi_S) d\xi_S. $$

By developing the likelihood $p(\mathbf{y} | \mathbf{X}_S, \xi_S)$ and the prior $p(\xi_S)$ in (13), we obtain

$$ \begin{aligned} p(\xi_S | \mathbf{X}_S, \mathbf{y}) & \propto \exp \left( -\frac{1}{2\sigma^2} [\mathbf{y} - \mathbf{X}_S \xi_S]^\top [\mathbf{y} - \mathbf{X}_S \xi_S] \right) \exp \left( -\frac{1}{2\tau^2} \xi_S^\top \mathbf{K}_S^{-1} \xi_S \right) \\ & \propto \exp \left( -\frac{1}{2} [\xi_S - \boldsymbol{\mu}_S]^\top \left[ \mathbf{X}_S^\top \mathbf{X}_S / \sigma^2 + \mathbf{K}_S^{-1} / \tau^2 \right] [\xi_S - \boldsymbol{\mu}_S] \right). \end{aligned} $$
Consequently, the posterior distribution {$\xi_S | X_S, y$} ~ $N(\mu_S, \Sigma_S)$, where

$$ \left\{ \begin{array}{l} \mu_S = \Sigma_S X_S^\top y / \sigma^2; \\ \Sigma_S^{-1} = X_S^\top X_S / \sigma^2 + K_S^{-1} / \tau^2. \end{array} \right. \qquad (14) $$

The formulas in (14), combined with the methodology developed in Section 4 and the smooth circulant embedding technique developed in Subsection 4.2, are highly useful for handling large datasets in multidimensional input spaces (see, Figure 8 in Example 5).

**Remark 1** (Improved Computational Efficiency in Posterior Inference).

* The matrix $\Sigma_S^{-1}$ in (14) can be obtained by applying the matrix inversion lemma (also known as the Woodbury formula) to $\tilde{\Sigma}_S$ in (12) (Williams and Rasmussen, 2006, Appendix A3):

$$ (T + UVW^\top)^{-1} = T^{-1} - T^{-1}U (V^{-1} + W^\top T^{-1}U)^{-1} W^\top T^{-1}, $$

assuming that all inverses exist.

* The predictive equations in (14) are significantly more efficient than the equivalent formulas in (12), presented in Lemma 6, especially for large datasets in multidimensional settings. Using the posterior precision matrix $\Sigma_S^{-1}$ is more advantageous than using $\tilde{\Sigma}_S^{-1}$, as it requires inverting only a $\prod_{\zeta=1}^d N_\zeta \times \prod_{\zeta=1}^d N_\zeta$ matrix, rather than both an $n \times n$ and a $\prod_{\zeta=1}^d N_\zeta \times \prod_{\zeta=1}^d N_\zeta$ matrix. It is worth noting that, when sampling from the posterior distribution (see Subsection 4.2), it is not necessary to compute the posterior mean $\mu_S$ or the posterior covariance matrix $\Sigma_S$ in (14). As a result, we only need to invert a single $\prod_{\zeta=1}^d N_\zeta \times \prod_{\zeta=1}^d N_\zeta$ matrix.

Now, let $\mathcal{C}$ be the set of functions that verify shape constraints such as monotonicity, convexity and boundedness, or any combination of these. The advantage of the model in (8) is that various shape constraints can be equivalently reformulated as linear inequality constraints on the Gaussian weight $\xi_S$, i.e.,

$$ \tilde{Y}_S \in \mathcal{C} \iff \xi_S \in \mathcal{E}(\xi_S), \qquad (15) $$

where $\mathcal{E}(\xi_S) = \{\xi_S \text{ s.t. } A\xi_S + b \ge 0\}$ is a polyhedron, with $A$ being the $m \times \prod_{\zeta=1}^d N_\zeta$ matrix of inequality constraints, and $b$ an $m$-dimensional vector. The matrix and vector of constraints, $A$ and $b$, depend on the shape constraints. For example, in one-dimensional case, if $\mathcal{C}$ is the set of monotone (non-decreasing) function, denoted by $C_m := \{f \in C^0(\mathcal{X}, \mathbb{R}) \text{ s.t. } f(x) \le f(y), \forall x < y\}$,
then $\tilde{Y}_N(\cdot)$, as defined in (5), belongs to $C_m$ if and only if $\xi \in \mathcal{E}(\xi)$, where

$$ \mathcal{E}(\xi) = \{\xi \in \mathbb{R}^N \text{ s.t. } \xi_j \le \xi_{j+1}, \forall j = 1, \dots, N-1\}. $$

In this case, we have

$$ A_{i,j} = \begin{cases} -1 & \text{if } j = i \quad \text{for any } i = 1, \dots, N-1; \\ 1 & \text{if } j = i+1 \quad \text{for any } i = 1, \dots, N-1; \\ 0 & \text{otherwise;} \end{cases} $$

and $b = [0, \dots, 0]^T \in \mathbb{R}^{N-1}$.

Before presenting the following result concerning monotonicity in the multidimensional case, let us recall that a function $f \in C([0, 1]^d, \mathbb{R})$ is said to be monotone (non-decreasing) with respect to all of its variables if $f(x) \le f(x')$, for any $x \le x' \in [0, 1]^d$, where $x \le x'$ denotes componentwise inequality, i.e., $x_1 \le x'_1, \dots, x_d \le x'_d$, with $x = [x_1, \dots, x_d]^T$.

**Proposition 1** (Monotonicity with respect to a subset of $\mathcal{X}$). If the finite-dimensional GP approximation is defined as in (8), then $\tilde{Y}_S(x_1, \dots, x_d)$ is monotone (non-decreasing) with respect to $x_\zeta$ if and only if

$$ Y(t_{j_1}^1, \dots, t_{j_{\zeta-1}}^{\zeta-1}, t_{j_\zeta}^{\zeta}, t_{j_{\zeta+1}}^{\zeta+1}, \dots, t_{j_d}^d) \le Y(t_{j_1}^1, \dots, t_{j_{\zeta-1}}^{\zeta-1}, t_{j_{\zeta+1}}^{\zeta}, t_{j_{\zeta+1}}^{\zeta+1}, \dots, t_{j_d}^d), $$

for any $j_\zeta = 1, \dots, N_\zeta - 1$ and for any $j_{\zeta^*} = 1, \dots, N_{\zeta^*}$, with $\zeta^* \in \mathcal{J}_{-\zeta} = \{1, \dots, d\} \setminus \{\zeta\}$. This leads to $m = (N_\zeta - 1) \prod_{\zeta^* \in \mathcal{J}_{-\zeta}} N_{\zeta^*}$ linear inequality constraints on the Gaussian weight $\xi_S$. In a similar way, we can express the monotonicity constraints on $\tilde{Y}_S$ with respect to any subset of $\mathcal{X}$. Finally, the monotonicity of $\tilde{Y}_S$ with respect to all of its input variables leads to $m = \sum_{\zeta=1}^d (N_\zeta - 1) \prod_{\zeta^* \in \mathcal{J}_{-\zeta}} N_{\zeta^*}$ linear inequality constraints.

As in the one-dimensional case, the posterior distribution $\{\tilde{Y}_S(\cdot) | y, \tilde{Y}_S \in \mathcal{C}\}$ can be derived from the tMVN distribution $\{\xi_S | X_S, y, \xi_S \in \mathcal{E}(\xi_S)\}$. Specifically, we have

$$ \{\xi_S | X_S, y, \xi_S \in \mathcal{E}(\xi_S)\} \sim \mathcal{T}\mathcal{N}(\mu_S, \Sigma_S; \mathcal{E}(\xi_S)) = \mathcal{T}\mathcal{N}(\mu_S, \Sigma_S; A, b), \quad (16) $$

where $\mathcal{T}\mathcal{N}(\mu_S, \Sigma_S; A, b)$ represents the tMVN distribution with mean vector $\mu_S$ and covariance matrix $\Sigma_S$ (as given in (14)), constraints matrix $A$, and constraints vector $b$. Thus, sampling
the finite-dimensional GP approximation (8) under linear inequality constraints necessitates sampling from a high-dimensional tMVN distribution, for which existing algorithms include Gibbs sampling (Taylor and Benjamini, 2016), Metropolis-Hastings (MH) (Cotter et al, 2013; Neal, 1999), the minimax tilting method accept-reject sampler (Botev, 2017), the Vecchia approximation of MVN (Cao and Katzfuss, 2023), and the highly efficient Hamiltonian Monte Carlo (HMC) sampler (Pakman and Paninski, 2014). If a stationary covariance function, such as the Matérn family (2), is employed, the uniform discretization of the input set $\mathcal{X}$ preserves the stationary property of the prior covariance matrix in (6). In this case, the prior covariance matrix $\tau^2 K_S$ exhibits Toepliz or block-Toeplitz structure, and efficient algorithms for generating the prior Gaussian vector $\xi_S$ are available, such as the FFT method developed by Wood and Chan (1994) and the fast large-scale algorithms proposed by Maatouk et al (2023, 2025c).

Before introducing the adaptive low-rank methodology and the efficient smooth circulant embedding technique for generating the posterior tMVN distribution (16), we define the two estimators used in this paper: mean a posteriori (mAP) and maximum a posteriori (MAP). The MAP estimator plays a crucial role in the low-rank methodology, the ordering of input variables, and the selection of hyperparameters. Unlike the mAP estimator, the MAP estimator is independent of the sampling procedure and is based solely on a quadratic optimization problem. Notably, the MAP estimator serves as a bridge between Bayesian inference and splines (Kimeldorf and Wahba, 1970). Recently, it has been utilized to address the problem of conditioning GPs with nonlinear datasets (Chen et al, 2025).

### Mean and maximum a posteriori (mAP and MAP)

In this paragraph, we define two posterior estimators based on the approach defined in (8): the mean a posteriori (mAP) and the maximum a posteriori (MAP). Their performance, in terms of prediction accuracy, is compared throughout the paper.

**Definition 1** (mAP estimator). The mean a posteriori (mAP) estimate of the Bayesian linear model $\tilde{Y}_S$ defined in (8), conditioned on linear inequality constraints (shape constraints) and noisy observations, is defined as follows:

$$ \bar{Y}_{\text{mAP}}^{S}(\boldsymbol{x}) := \text{E} \left[ \tilde{Y}_{S}(\boldsymbol{x}) | \boldsymbol{y}, \tilde{Y}_{S} \in \mathcal{C} \right] = \Phi_{S}(\boldsymbol{x})^{\top} \bar{\mu}_{S}, \quad \boldsymbol{x} \in \mathcal{X}, $$
where $\bar{\mu}_S := E[\xi_S|y, \xi_S \in \mathcal{E}(\xi_S)]$ is the posterior mean which is computed from Monte Carlo simulations, with $\mathcal{E}(\xi_S)$ the set of linear inequality constraints on the $\prod_{\zeta=1}^d N_\zeta$-dimensional Gaussian weights vector $\xi_S$.

**Definition 2 (MAP estimator).** The maximum a posteriori (MAP) estimate of the Bayesian linear model $\tilde{Y}_S$ defined in (8), conditioned on linear inequality constraints (shape constraints) and noisy observations, is defined as follows:

$$ \hat{Y}_{\text{MAP}}^{S}(\boldsymbol{x}) := \Phi_{S}(\boldsymbol{x})^{\top} \hat{\mu}_{S}, \quad \boldsymbol{x} \in \mathcal{X}, $$

where $\hat{\mu}_S$ is the posterior mode computed by the quadratic problem (QP), see (18) below.

Let $\hat{\mu}_S$ be the posterior mode that maximises the probability density function (pdf) of the posterior tMVN distribution in (16). Then, maximizing this posterior pdf is equivalent to maximize the following QP

$$ \hat{\mu}_S := \arg\max_{\xi_S \text{ s.t. } A\xi_S + b \ge 0} \{-[\xi_S - \mu_S]^\top \Sigma_S^{-1} [\xi_S - \mu_S]\}, \quad (17) $$

with condition mean $\mu_S$ and conditional covariance $\Sigma_S$ as in (14). This is a QP subject to linear inequality constraints (Boyd and Vandenberghe, 2004; Goldfarb and Idnani, 1983). The optimization problem in (17) is equivalent to

$$ \hat{\mu}_S := \arg\min_{\xi_S \text{ s.t. } A\xi_S + b \ge 0} \left\{ \xi_S^\top \Sigma_S^{-1} \xi_S - 2y^\top X_S \xi_S / \sigma^2 \right\}. \quad (18) $$

It is worth noting that the posterior mode $\hat{\mu}_S$ is used as a starting point for the proposed MCMC sampler developed in Subsection 4.2. Furthermore, it plays an important role in avoiding numerical instability in the sampling procedure for a large value of the approximation parameter $\eta$ (see Maatouk et al (2024) and Subsection 4.2 below for more details). Finally, it is an efficient estimate point as it is based on a QP and on the formulas in (14).

## Comments on the MAP

* The MAP estimator $\hat{Y}_{\text{MAP}}^S(\cdot)$ is independent of the sampling procedure and is determined by solving a QP subject to linear inequality constraints, as given in (18). However, it does depend on the noise variance parameter $\sigma^2$ and the discretization $S$ of the input set $\mathcal{X}$, since both $\Phi_S$ and $\hat{\mu}_S$ are functions of the discretization.
* The asymptotic behavior of this estimator has been discussed in both Bay et al (2016) and Grammont et al (2024), where the authors established a generalization of the Kimeldorf-Wahba correspondence (Kimeldorf and Wahba, 1970) for both constrained noisy and noise-free cases. This generalization extends the original correspondence between Bayesian estimation and optimal interpolation in a Reproducing Kernel Hilbert Space (RKHS) to include linear inequality constraints such as boundedness, monotonicity, or convexity.

* In the present paper, the MAP estimator $\hat{Y}_{\text{MAP}}^S(\cdot)$ plays an important role in the low-rank methodology and the ordering of input variables (Section 4 below). This is due to the fact that its computational complexity is much lower than that of the mAP estimator $\overline{Y}_{\text{mAP}}^S(\cdot)$.

* The efficient formulas (14) are used to compute the posterior mode (18) and, consequently, the MAP estimator in Definition 2, thereby reducing the computational complexity to $O((\prod_{\zeta=1}^d N_\zeta)^3)$ instead of $O((\prod_{\zeta=1}^d N_\zeta)^3 + n^3)$. It is worth noting that, when computing the posterior mode (18), it is not necessary to compute the posterior mean $\mu_S$ or the posterior covariance matrix $\Sigma_S$ in (14). As a result, we only need to invert a single $\prod_{\zeta=1}^d N_\zeta \times \prod_{\zeta=1}^d N_\zeta$ matrix.

As explained before (16), in higher dimensions (i.e., when $\mathcal{X} = [0, 1]^d$ and $d \gg 2$), employing a uniform discretization of the input set $\mathcal{X}$ requires sampling a truncated Gaussian vector of dimension $\prod_{\zeta=1}^d N_\zeta$, where $N_\zeta$ denotes the number of knots for the input dimension $\zeta \in \{1, \dots, d\}$. When $N_\zeta$ and $d$ are large, this approach becomes computationally prohibitive and impractical. Conversely, in low dimensions with uniform subdivision, for any $d \ge 1$, when the values of $\{N_\zeta\}$, for $\zeta \in \{1, \dots, d\}$ are too small, the model in (8) yields poor predictions, as illustrated in the top-left panel of Figure 1. To address these issues, we develop a method that combines low-rank with input variable ordering. This approach is based on iteratively adding optimal, non-uniform, one-dimensional knots and performing the ordering of active input variables. It is well-suited to multidimensional input settings and can efficiently handle a large number of observations. This efficiency stems from the fact that the developed MCMC approach based on Elliptical Slice Sampling (ESS), see Subsection 4.2, avoids matrix inversion and does not require computing the full posterior distribution (12) or (14) during the sampling procedure of the tMVN in (16).
# 4 Low-rank methodology and variable ordering

## 4.1 Adaptive discretization methodology

As shown in Definitions 1 and 2, the computational cost of the two proposed estimators mAP and MAP (denoted by $\bar{Y}_{\text{mAP}}^S$ and $\hat{Y}_{\text{MAP}}^S$, respectively), depends on the discretization of the input set $\mathcal{X}$. In this section, we develop a low-rank strategy and the ordering of active input variables for the finite-dimensional GP approximation (8). Low-rank refers here to reducing the dimensionality of the Gaussian vector $\xi_S$ in (8) by computing an optimal non-uniform discretization of $\mathcal{X}$. The proposed approach is applicable to both low- and high-dimensional input settings and can handle both small and large datasets. The proposed methodology is based on the iterative construction of optimal non-uniform one-dimensional knots and the ordering of input variables. Unlike the mAP estimator, the MAP estimator is independent of the sampling procedure and only depends on a quadratic optimization problem (18). For this reason, the MAP estimator plays a key role in our methodology. Additionally, we use the predictive equations in (14), as they are more efficient for handling large datasets than those in (12) (i.e., when $n \gg \prod_{\zeta=1}^d N_\zeta$). For simplicity of presentation, we denote by $\zeta$ the variable representing the index of an input dimension. Consequently, $\zeta \in \{1, \dots, d\}$, where $d \ge 1$ denotes the number of input dimensions.

To describe the procedure of the proposed algorithm, let us first recall that $\mathcal{S} = (s^{(1)}, \dots, s^{(d)})$ denotes the subdivision of the input set $\mathcal{X}$. For example, for any $\zeta \in \{1, \dots, d\}$, we have $s^{(\zeta)} = \{t_1^{\zeta}, \dots, t_{N_\zeta}^{\zeta}\}$, representing a given ordered partition of the input dimension $\zeta$ at a given iteration. The aim is to enrich $\mathcal{S}$ at each step of the proposed sequential procedure, starting from a basic uniform partition (see the sequential procedure below). It is worth noting that this basic uniform partition depends on the shape constraints (Remark 2). Hereafter, we develop the sequential procedure:

* **Phase 0 - Initialization:** The first step of the proposed approach involves constructing a square uniform discretization of the input set $\mathcal{X}$. Specifically, we fix $N_\zeta^* = 2$, for any $\zeta \in \{1, \dots, d\}$. Consequently, for any $\zeta$, the discretization set is given by $s_\star^{(\zeta)} = s^* = \{t_1^\zeta, t_{N_\zeta}^\zeta\} = \{0, 1\}$. As a result, the initial subdivision of $\mathcal{X}$ is defined as $\mathcal{S}^* = \{s^*, \dots, s^*\}$. In this particular case, the covariance matrix of the basis coefficients, $\tau^2 K_{\mathcal{S}^*}$, is block-Toeplitz. To compute optimal one-dimensional knots, we define a knot set $s_{\text{fix}}$, from which the optimal knots are selected. The goal is to sequentially compute the optimal non-uniform one-dimensional knots chosen from $s_{\text{fix}}$. This one-dimensional uniform discretization $s_{\text{fix}}$ can be made as fine as desired. However, based on several testing scenarios, excessively
fine subdivisions tend to produce a biased model. To mitigate this effect, we impose a minimum distance of 0.05 between any two optimal one-dimensional knots throughout this work. In other words, newly inserted knots must not be too close to existing ones.

* **Phase 1 - Active variable ordering:** In this step, we determine the ordering of the input variables $\zeta \in \{1, \dots, d\}$. Specifically, we select the first active input variable $\zeta \in \{1, \dots, d\}$ such that adding an optimal one-dimensional knot in this direction optimizes a given criterion of the proposed model. We continue selecting active input variables in this manner until a fixed optimal ordering is obtained. This strategy is advantageous because the active input variables remain ordered throughout all iterations, which significantly reduces the computational complexity of the proposed algorithm. In the present paper and without loss of generality, we decide using the mean squared error (MSE) criterion through the MAP estimator $\hat{Y}_{MAP}^S$ (Definition 2):

$$ \text{MSE}_{\text{MAP}}(S) = \frac{1}{n} \sum_{i=1}^{n} \left( \hat{Y}_{\text{MAP}}^{S}(\mathbf{x}_i) - y_i \right)^2 . $$

The MSEMAP(S) is an interesting criterion, as it depends on the observed data $\{y_i\}$ and on the efficient MAP estimator. The MAP estimator $\hat{Y}_{MAP}^S$ is independent of the sampling procedure, and computing it only requires solving the quadratic optimization problem in (18). Clearly, it depends on the subdivision S of $\mathcal{X}$, as the posterior mode $\hat{\mu}_S$ in (18) itself depends on S.

* **Phase 2 - Refining input direction $\zeta$:** For the current input direction $\zeta$ of the optimal ordering, we denote by $s_{old}^{(\zeta)}$ a non-uniform discretization consisting of a set of optimal one-dimensional knots associated with the input direction $\zeta \in \{1, \dots, d\}$. Refinement in the direction $\zeta$ is performed by adding a new knot $t \in \{s_{fix} \setminus s_{old}^{(\zeta)}\} \subset [0, 1]$. The resulting refined discretization $s_{new}^{(\zeta)}$ is then defined as:

$$ \left\{ \begin{array}{l} s_{new}^{(\zeta)} = \text{ord}(s_{old}^{(\zeta)} \cup t); \\ N_{\zeta}^{new} = N_{\zeta}^{old} + 1; \end{array} \right. $$

where $N_{\zeta}^{old}$ represents the number of non-uniform one-dimensional knots at a fixed iteration for a given input variable $\zeta \in \{1, \dots, d\}$. More precisely, the new one-dimensional knot $t \in \{s_{fix} \setminus s_{old}^{(\zeta)}\}$ is selected by minimizing the MSEMAP(S) criterion. Another way is to formulate an optimization problem to determine the optimal non-uniform one-dimensional knots, which can again be achieved by using the MAP estimator as the predictive function.
Thus,

$$S_{\text{new}} = (s^{(1)}, \dots, s^{(\zeta-1)}, \underbrace{\text{ord}(s_{\text{old}}^{(\zeta)} \cup t), s^{(\zeta+1)}, \dots, s^{(d)})_{s_{\text{new}}^{(\zeta)}}, \dots, s^{(\zeta-1)}, \dots, s^{(1)}) \quad (19)$$

where ord(·) is function that ordering the one-dimensional knots in an increasing order. For instance, if $t \in [0, 1]$ is such that $t_q^\zeta < t < t_{q+1}^\zeta$ and $t_q^\zeta < t_{q+1}^\zeta$ for any $q = 1, \dots, N_\zeta - 1$, then, $\text{ord}(\{t_{j_\zeta}^\zeta\} \cup t) = (t_1^\zeta, \dots, t_q^\zeta, t, t_{q+1}^\zeta, \dots, t_{N_\zeta}^\zeta)$, where $\{t_{j_\zeta}^\zeta\}$ represents $s^{(\zeta)}$ in (19). We continue adding one-dimensional knots in this direction and updating the corresponding one-dimensional subdivision $s_{\text{new}}^\zeta$ until either the stopping criterion of the proposed algorithm is met or the maximum number of one-dimensional knots, $N_\zeta$, is reached. It is worth noting that, in all numerical examples presented in this paper, the option to impose a maximum number of one-dimensional knots is not employed.

It is worth noting that during the three aforementioned steps (initialization, the ordering of active input variables, and input direction refinement), the constraint matrix $A$ and vector $b$, which define the set of linear inequality constraints $\mathcal{E}(\xi_S)$ on the basis coefficients $\xi_S$, as well as the design matrix $X_S$ and the prior covariance matrix $K_S$, must be updated at each iteration. As previously mentioned, the proposed approach avoids computing the full posterior distribution given in (14) during the sampling procedure. Instead, it requires only the computation of the matrix $\Sigma_S^{-1} = X_S^\top X_S + K_S^{-1}$ (needed for the posterior mode $\hat{\mu}_S$ in (18)) at each iteration of the procedure. The proposed sequential algorithm minimizes a given criterion by selecting either an input variable to activate or a one-dimensional knot to insert. In this work, the chosen criterion is the MSE$_{MAP}(S)$, which is based on the MAP estimator $\hat{Y}_{MAP}^S$ and the given data. The proposed algorithm is outlined in Algorithm 1. For clarity of presentation, certain practical optimizations, such as avoiding redundant computations, are not shown.

**Remark 2.** If convexity (or concavity) constraints are required, then, the subdivision initialization $s^*$ should be updated. In the algorithm developed in this paper, we suggest starting with $s^* = \{0, 0.5, 1\}$ and refining the corresponding input set if needed. For simplicity, let us consider a one-dimensional input space as an example. In this case, the input set $X$ is the unit interval containing $N \ge 3$ non-uniform knots $\{t_j\}$. Then, the finite-dimensional GP approximation $\tilde{Y}_N$ is convex in the entire domain if and only if $\xi \in \mathcal{E}(\xi)$, where

$$\mathcal{E}(\xi) = \left\{ \xi \in \mathbb{R}^N : \frac{\xi_j - \xi_{j-1}}{t_j - t_{j-1}} \le \frac{\xi_{j+1} - \xi_j}{t_{j+1} - t_j}, j = 2, \dots, N-1 \right\},$$
**Algorithm 1** Low-rank and input variable ordering for large datasets

Initialization: fix and compute the following terms:

* Input: $N_ζ$, $S^*$, $y$, and $s_{fix}$

* Compute the terms: $X_{S^*}$, $\mathcal{E}(\xi_{S^*})$, and $\text{MSE}_{\text{MAP}}(S^*)$

Phase 1 – Active variable ordering

Dim $\leftarrow \{1, \dots, d\}$

for Ord from 1 to d do

Find $\zeta$ in Dim, which minimizes $\min_{t \in s_{fix}} \text{MSE}_{\text{MAP}}(S^* \cup t)$

Dim $\leftarrow \text{Dim} \setminus \{\zeta\}$

VarOrd[Ord] $\leftarrow \zeta$

end for

Phase 2 – Refining input direction

for each $\zeta$ in the optimal ordering variables VarOrd do

while $\text{MSE}_{\text{MAP}}(S^*)$ is sufficiently improved do

Add a one-dimensional knot $t \in \{s_{fix} \setminus s^*\}$ that minimizes $\text{MSE}_{\text{MAP}}(S^* \cup t)$:

$$S^* \leftarrow S^* \cup t := (s^*_{(1)}, \dots, s^*_{(\zeta-1)}, \text{ord}(s^*_{(\zeta)} \cup t), s^*_{(\zeta+1)}, \dots, s^*_{(d)})$$

$$N_ζ \leftarrow N_ζ + 1$$

Update $X_{S^*}$, $\mathcal{E}(\xi_{S^*})$, and $\text{MSE}_{\text{MAP}}(S^*)$

end while

end for

Phase 3 – Sampling the posterior distribution

Sample the posterior tMVN $\{\xi_{S^*} | X_{S^*}, y, \xi_{S^*} \in \mathcal{E}(\xi_{S^*})\}$ using the MCMC-ESS approach developed in Subsection 4.2.

with $\xi = [\xi_1, \dots, \xi_N]^\top$. Finally, let us mention that this initialization $s^* = \{0, 0.5, 1\}$ should be taken into account when convexity or concavity constraints are among multiple shape constraints.

## 4.2 Posterior distribution using Elliptical Slice Sampling

In this section, we extend the efficient smooth circulant embedding technique proposed by Ray et al (2020) to support implementation of the proposed approach in multidimensional input spaces. This technique is based on a smooth relaxation of the set of constraints and on the use of the ESS (Murray et al, 2010). We employed this approach because it demonstrates excellent flexibility and high accuracy in the context of Bayesian shape-restricted function estimation (Maatouk et al, 2025b, 2024, 2025c). We explore its application to the adaptive low-rank methodology in both low- and high-dimensional input spaces, as well as the ordering of input variables. Furthermore, we discuss its capability to address hyperparameter estimation in these cases.
According to both (15) and (16), we are interested in sampling from a constrained posterior distribution with the following form:

$$p(\xi_S | \mathbf{X}_S, \mathbf{y}, \xi \in \mathcal{E}(\xi_S)) \propto \underbrace{\exp\left(-\frac{1}{2\sigma^2}\|\mathbf{y} - \mathbf{X}_S\xi_S\|^2\right)}_{\text{likelihood function}} \underbrace{\exp\left(-\frac{1}{2\tau^2}\xi_S^\top \mathbf{K}_S^{-1} \xi_S\right)}_{\text{(untruncated) Gaussian prior}} \underbrace{\mathbf{1}_{\mathcal{E}(\xi_S)}(\xi_S)}_{\text{constraint}}, \quad (20)$$

where $\mathbf{y} \in \mathbb{R}^n$ is the vector of data, $\mathbf{X}_S$ is the $n \times \prod_{\zeta=1}^d N_\zeta$ design matrix, $\mathcal{E}(\xi_S)$ is a convex set representing the set of linear inequality constraints on the $\prod_{\zeta=1}^d N_\zeta$-dimensional Gaussian weight $\xi_S$ with $\mathbf{K}_S$ a positive-definite matrix. In Maatouk et al (2025b, 2024, 2025c); Ray et al (2020), the smooth circulant embedding technique developed in this section has been employed in cases where the prior covariance matrix $\tau^2 \mathbf{K}_S$ is Toeplitz. In the present paper, this property no longer holds because the discretization of the input set $\mathcal{X}$ is no longer uniform. The main idea is to approximate the indicator function in (20) by a smooth function and incorporate it into the likelihood function. Let us recall the following result proved by Maatouk et al (2024), which generalizes the results in Ray et al (2020):

**Proposition 2** (Approximation of indicator function of linear inequality constraints). Let $\mathbf{a}$ be a $\prod_{\zeta=1}^d N_\zeta$-dimensional vector and $b \in \mathbb{R}$. Then, the indicator function $\mathbf{1}_{\{\mathbf{a}^\top \xi_S + b \ge 0\}}$ can be approximated by a scaled logistic sigmoid function as follows:

$$\mathbf{1}_{\{\mathbf{a}^\top \xi_S + b \ge 0\}} \approx \frac{1}{1 + \exp[-\eta(\mathbf{a}^\top \xi_S + b)]}, \quad (21)$$

where the parameter $\eta > 0$ controls the quality of the approximation. As $\eta$ increases, the quality of the approximation improves. Furthermore, we have

$$\left| \mathbf{1}_{\{\mathbf{a}^\top \xi_S + b \ge 0\}} - \frac{1}{1 + e^{-\eta(\mathbf{a}^\top \xi_S + b)}} \right| \le \frac{1}{1 + e^{\eta|\mathbf{a}^\top \xi_S + b|}},$$

for any $\prod_{\zeta=1}^d N_\zeta$-dimensional vector $\xi_S$.

By substituting the proposed approximation (21) in (20), we get:

$$\tilde{p}_\eta(\xi_S | \mathbf{X}_S, \mathbf{y}, \xi_S \in \mathcal{E}(\xi_S)) \propto \underbrace{\exp\left(-\frac{1}{2\sigma^2}\|\mathbf{y} - \mathbf{X}_S\xi_S\|^2\right)}_{\text{new likelihood function}} \times \prod_{k=1}^m \frac{1}{1 + e^{-\eta(a_k^\top \xi_S + b_k)}} \\ \times \underbrace{\exp\left(-\frac{1}{2\tau^2}\xi_S^\top \mathbf{K}_S^{-1} \xi_S\right)}_{\text{Gaussian prior}}, \quad (22)$$
where $\boldsymbol{\alpha}_k$ is a $\prod_{\zeta=1}^d N_\zeta$-dimensional vector representing the $k^{\text{th}}$ row of the constraint matrix $\mathbf{A}$ and $b_k$ is the $k^{\text{th}}$ component of the vector $\boldsymbol{b} \in \mathbb{R}^m$. The above approximate posterior pdf (22) is proportional to a product of a new likelihood function and a prior MVN distribution:

$$ \tilde{p}_{\eta}(\boldsymbol{\xi}_S | \mathbf{X}_S, \mathbf{y}, \boldsymbol{\xi}_S \in \mathcal{E}(\boldsymbol{\xi}_S)) \propto L(\cdot) \mathcal{N}(\boldsymbol{\xi}_S; \mathbf{0}, \tau^2 \mathbf{K}_S). \quad (23) $$

In this context, the ESS can be employed which is seen as a parameterization of Metropolis-Hasting (Neal, 1999):

$$ \boldsymbol{\xi}'_S = \sin(\theta) \boldsymbol{\nu}_S + \cos(\theta) \boldsymbol{\xi}_S, \quad \boldsymbol{\nu}_S \sim \mathcal{N}(\mathbf{0}, \tau^2 \mathbf{K}_S). $$

where the angle $\theta$ is uniformly generated from a $[\theta_{\min}, \theta_{\max}]$ interval, which shrinks exponentially fast until an acceptable state is reached (see Murray et al (2010) for more detailed guidelines). This is a variant of the so-called preconditioned Crank-Nicolson (pCN) method (Cotter et al, 2013). The pCN algorithm mixes the current state with a Gaussian sample drawn from the prior, leading to dimension-independent acceptance rates and robust performance in Bayesian inference problems with Gaussian priors. The logarithm of the new likelihood in (22) is evaluated at each MCMC iteration and it is expressed as follows:

$$ \log[L(\boldsymbol{\xi}_S)] = -\frac{1}{2\sigma^2} \|\mathbf{y} - \mathbf{X}_S \boldsymbol{\xi}_S\|^2 - \sum_{k=1}^{m} \log \left(1 + e^{-\eta(\boldsymbol{\alpha}_k^\top \boldsymbol{\xi}_S + b_k)}\right). \quad (24) $$

The computational complexity for evaluating the log-likelihood in (24) at each MCMC iteration is $O(\max(n \prod_{\zeta=1}^d N_\zeta, m \prod_{\zeta=1}^d N_\zeta))$, where $n$ is the number of samples, $m$ is the number of constraints, and $\prod_{\zeta=1}^d N_\zeta$ denotes the dimension of the MVN prior (i.e., the number of knot points).

**Proposition 3** (Boundedness constraints in multidimensional cases). If $\mathcal{C}$ is the set of functions verifying boundedness constraints

$$ \mathcal{C} = \{f \in C^0(\mathcal{X}, \mathbb{R}) \text{ s.t. } l \le f(\boldsymbol{x}) \le u, \boldsymbol{x} \in \mathcal{X}, l < u \in \mathbb{R}\}, \quad (25) $$

then, according to the notation in (15), the corresponding set of linear inequality constraints on the basis coefficients is $\mathcal{E}(\boldsymbol{\xi}_S) = \{\boldsymbol{\xi}_S \text{ s.t. } l \le \boldsymbol{\xi}_S \le \mathbf{u}\}$, where $l = [l, \dots, l]^\top$ and $\mathbf{u} = [u, \dots, u]^\top$ are the $\prod_{\zeta=1}^d N_\zeta$-dimensional lower and upper bound vectors, respectively. Moreover, the matrix
and vector of linear inequality constraints are $\mathbf{A} = [I_{\prod_{\zeta=1}^d N_\zeta}, -I_{\prod_{\zeta=1}^d N_\zeta}]^\top$ and $\mathbf{b} = [-\mathbf{l}^\top, \mathbf{u}^\top]^\top$, of dimension $2\prod_{\zeta=1}^d N_\zeta \times \prod_{\zeta=1}^d N_\zeta$ and $2\prod_{\zeta=1}^d N_\zeta$, respectively. Let us recall that $I_{\prod_{\zeta=1}^d N_\zeta}$ is the $\prod_{\zeta=1}^d N_\zeta \times \prod_{\zeta=1}^d N_\zeta$ identity matrix. This involves $m = 2\prod_{\zeta=1}^d N_\zeta$ linear inequality constraints on the basis coefficients $\mathbf{A}\boldsymbol{\xi}_S + \mathbf{b} \ge 0$. For simplicity, the boundedness constraints for the single-index $\zeta \in \{1, \dots, d\}$ involves a linear inequality constraints with a $2N_\zeta \times N_\zeta$ matrix $[I_{N_\zeta}, -I_{N_\zeta}]^\top$ and a $2N_\zeta$-dimensional vector $[-l, \dots, -l, u, \dots, u]^\top$.

**Corollary 2 (Generalization of boundedness constraints).** In line with Proposition 3, let $\mathcal{C}$ be the set of functions defined as

$$ \mathcal{C} = \{ f \in C^0(\mathcal{X}, \mathbb{R}) \text{ s.t. } f_l(\mathbf{x}) \le f(\mathbf{x}) \le f_u(\mathbf{x}), \mathbf{x} \in \mathcal{X} \}, $$

where $f_l$ and $f_u$ are continuous lower and upper bounding functions on $\mathcal{X}$, respectively, that define a convex set. Then,

$$ \tilde{Y}_S \in \mathcal{C} \iff f_l(t_{j_1}^1, \dots, t_{j_d}^d) \le \xi_{j_1, \dots, j_d} \le f_u(t_{j_1}^1, \dots, t_{j_d}^d), $$

for any $j_\zeta = 1, \dots, N_\zeta$, with $\zeta \in \{1, \dots, d\}$ (see Maatouk et al (2025a) for more details). In fact, $f_l(t_{j_1}^1, \dots, t_{j_d}^d)$ and $f_u(t_{j_1}^1, \dots, t_{j_d}^d)$ play the role of the lower and upper bound vectors $\mathbf{l}$ and $\mathbf{u}$, respectively, in Proposition 3.

**Lemma 7 (A simple likelihood formulation for boundedness constraints).** If the Gaussian vector $\boldsymbol{\xi}_S$ is bounded (i.e., belongs to $\mathcal{E}(\boldsymbol{\xi}_S)$) as in Proposition 3, then the sum in the log-likelihood function (24) has the following form:

$$ \sum_{k=1}^{2\prod_{\zeta=1}^{d} N_{\zeta}} \log \left( 1 + e^{-\eta(\boldsymbol{a}_k^{\top} \boldsymbol{\xi}_S + b_k)} \right) = \sum_{k=1}^{\prod_{\zeta=1}^{d} N_{\zeta}} \log \left( 1 + e^{-\eta(\boldsymbol{\xi}_S - \mathbf{l})_k} + e^{-\eta(\mathbf{u} - \boldsymbol{\xi}_S)_k} + e^{-\eta(\mathbf{u} - \mathbf{l})_k} \right). \quad (26) $$

It is worth noting that the expression in (26) eliminates the necessity of performing matrix products and hence, reducing the computational complexity to $\mathcal{O}(\prod_{\zeta=1}^d N_\zeta)$ instead of $\mathcal{O}((\prod_{\zeta=1}^d N_\zeta)^2)$. This significant advantage is explored when boundedness constraints are required for both individual and multiple shape constraints in both one-dimensional and multi-dimensional input settings (Examples 3 and 4). It is worth noting that, based on the result provided in Lemma 7, generating 20,000 samples from a tMVN distribution of dimension $d = 5,000$ is completed in three minutes (Maatouk et al, 2024).
## 4.3 Hyperparameters selection and full conditionals

Consider the Bayesian linear model in (8) and the prior information in Subsection 3.3. The noise variance and the signal variance hyperparameters, which appear in the regression model (1) and the prior covariance function (2), respectively, are generally unknown and must be estimated. In this context, the joint posterior distribution is given by:

$$ \pi(\mathbf{y}, \mathbf{X}_S, \boldsymbol{\xi}_S, \sigma^2, \tau^2) \propto \frac{1}{(\sigma^2)^{n/2}} \exp\{-\frac{1}{2\sigma^2} \|\mathbf{y} - \mathbf{X}_S\boldsymbol{\xi}_S\|_2^2\} \\ \times \frac{1}{(\tau^2)^{\frac{1}{2}} \prod_{\zeta=1}^d N_\zeta} \exp\{-\frac{1}{2\tau^2} \boldsymbol{\xi}_S^\top \mathbf{K}_S^{-1} \boldsymbol{\xi}_S\} \mathbf{1}_{\mathcal{E}(\boldsymbol{\xi}_S)}. $$

Then, {$\boldsymbol{\xi}_S | \mathbf{X}_S, \mathbf{y}, \sigma^2, \tau^2$} is a MVN distribution truncated on $\mathcal{E}(\boldsymbol{\xi}_S)$. Additionally, we have

$$ \{\sigma^2 | \mathbf{y}, \boldsymbol{\xi}_S, \tau^2\} \sim \mathcal{IG}(n/2, \|\mathbf{y} - \mathbf{X}_S \boldsymbol{\xi}_S\|_2^2/2); $$

$$ \{\tau^2 | \mathbf{y}, \boldsymbol{\xi}_S, \sigma^2\} \sim \mathcal{IG}(\frac{1}{2} \prod_{\zeta=1}^d N_\zeta, \boldsymbol{\xi}_S^\top \mathbf{K}_S^{-1} \boldsymbol{\xi}_S); $$

where $\mathcal{IG}$ represents the inverse Gamma distribution. Consequently, a Gibbs sampler is employed to draw samples from the joint posterior distribution of {$\boldsymbol{\xi}_S, \tau^2, \sigma^2$}.

# 5 Numerical examples

The aim of this section is to evaluate the performance of the proposed methodology using synthetic datasets and various shape constraints. It presents numerical examples in both one-dimensional and multidimensional cases, addressing individual, multiple, and complex shape constraints. Additionally, the capability of the approach to handle large datasets is analyzed. Finally, the performance of the mAP estimator (Definition 1) compared to the MAP estimator (Definition 2) in terms of prediction accuracy is examined. In all the numerical examples in this section, the shape constraints on the target function are treated as known prior information.

## 5.1 One-dimensional input space (complex shape constraints)

In this section, we illustrate the performance of the proposed methodology using two challenging examples in a one-dimensional input space, i.e., $d = 1$.

**Example 2 (Extreme case).** The aim of this section is to demonstrate the flexibility of the proposed approach in handling extreme cases. To do this, we consider the following target
function

$$f(x) = \begin{cases} x \sin(\frac{7}{3}\pi(x - 0.3)) & \text{if } x \in [0, 0.8]; \\ -0.4 & \text{otherwise.} \end{cases}$$

This situation is challenging because $f$, bounded between $-0.4$ and $0.6$, remains flat in the interval $[4/5, 1]$ and reaches the lower bound. To make the situation more challenging, we assume that the $n = 100$ training samples (black stars) are generated from (1) using the trigonometric function $x \sin(\frac{7}{3}\pi(x - 0.3))$ with a standard deviation of $\sigma = 0.05$. This implies that the observations (black stars in Figure 4) do not follow the target function's behavior in the second part of the domain, $[4/5, 1]$. Consequently, some training samples are considered outliers because they are not generated by the flat function and do not adhere to bounded constraints.

<figure><img src="image_6.png" /><figcaption>Fig. 4 Bounded function estimation in an extreme case.</figcaption></figure>

In Figure 4, we observe that the proposed approach splits the input domain $\mathcal{X}$ into $N = 8$ non-uniform knots (gray vertical dashed lines) at the locations where the function changes behavior. Additionally, we observe that the posterior distribution (gray shaded area, MAP and mAP) fits the data well in the first part of the domain $[0, 4/5]$, respects boundedness constraints across the entire domain, and effectively captures the flat behavior of the underlying function in the second part of the domain $[4/5, 1]$. This is due to the excellent flexibility of the proposed approach, whereas other MCMC samplers, such as HMC (Pakman and Paninski, 2014), fail to generate the posterior distribution in this extreme case. This may be because these approaches rigidly follow the behavior of the data, limiting their effectiveness.

**Example 3 (Multiple shape constraints).** In this numerical example, we show the effectiveness of the proposed methodology to handle cases when multiple shape constraints are required together. To do this, we consider the function $x \mapsto 1/\{1 + \exp(-10x)\}$. This function is monotone (non-decreasing) and bounded from above by 1. According to (15), the boundedness and monotonicity (non-decreasing) constraints on the Gaussian vector $\xi_S$ translate into the
following linear inequalities:

$$ \xi_j \le 1, \quad \forall j = 1, \dots, N \quad \text{and} \quad \xi_1 \le \dots \le \xi_N. $$

Together, these conditions result in only $m = N$ effective inequality constraints, namely $\xi_1 \le \dots \le \xi_N \le 1$.

<figure><img src="image_7.png" /><figcaption>Fig. 5 Performance illustration of constrained GP approximation under multiple shape constraints (monotonicity and boundedness), with an optimal non-uniform discretization (left panel) and a uniform discretization (right panel). The underlying function is $x \mapsto 1/(1 + \exp(-10x))$ which is monotone increasing on $[0, 1]$ and bounded from above by 1.</figcaption></figure>

In the left panel of Figure 5, the proposed methodology is employed, while in the right panel, a naive uniform discretization is used to place the $N = 50$ equally spaced knots. The developed methodology successfully handles multiple shape constraints and places optimal non-uniform knots to obtain a “good” approximation of the target function. It drastically reduces computational complexity, as it requires only a few non-uniform knots to accurately fit the given data (black stars). Furthermore, a large number of equally spaced knots leads to a poor posterior distribution (gray shaded area in the right panel), and thus, to poor predictive models. This is clearly related to the mass-shifting phenomenon, which is more pronounced in high dimensions (i.e., when the number of knots $N$ is large) and when the correlation between the components of the tMVN distribution is more significant. It is worth noting that, as mentioned previously, the MAP estimator is not affected by this phenomenon, as it is independent of the sampling procedure (i.e., sampling the posterior tMVN distribution (16)).

## 5.2 Multidimensional input space with large datasets

In this section, we investigate the feasibility of the proposed methodology in a multidimensional input space with a large dataset. Throughout this section, we use a tensor product of 1D Matérn covariance kernel (2) is used with a smoothness parameter $\nu = 5/2$.
**Example 4 (Boundedness in two-dimensional case).** We consider the target function

$$f(x_1, x_2) = \frac{1}{[1 + (10x_1)^4]} + \frac{1}{2} \exp[-100(x_1 - 0.5)^2] + x_2^2, \quad \forall (x_1, x_2) \in \mathcal{X},$$

which presents a challenging scenario due to its higher variability with respect to the first input variable. Consequently, we anticipate that the proposed algorithm will adaptively allocate more knots along the $x_1$-dimension. The target monotone function $f$ is evaluated on a Latin-hypercube design at 900 locations using `lhsDesign` function implemented in the R-package `DiceDesign`. The training data $\mathbf{y}$ are generated from (1) using the target function $f$ and a noise standard deviation of 0.01% of the range of $\mathbf{y}$.

<figure><img src="image_8.png" /><figcaption>**Fig. 6** Bounded function estimation for the two-dimensional case. Top: side-by-side comparison of the two estimators, mAP and MAP, with the true function. Bottom: the three top functions at $x_2 = 0$ (left) and $x_1 = 0$ (right). The gray vertical dashed lines represent the locations of the non-uniform one-dimensional knots.</figcaption></figure>

In the top panels of Figure 6, we present side-by-side comparisons of the two proposed estimators, mAP and MAP, along with the true monotone function. The proposed approach effectively captures the behavior of the target function and respect boundedness constraints in the entire domain. In the bottom panels, we illustrate the three top functions at $x_2 = 0$ (left) and $x_1 = 0$ (right). The gray vertical dashed lines indicate the locations of the non-uniform (one-dimensional) knots. Notably, the proposed approach prioritizes finding optimal non-uniform knots in the $x_1$-dimension, which exhibits greater variation. According to the notation introduction in Section 4, we obtain $N_1 = 13$ and $N_2 = 3$ optimal non-uniform (one-dimensional)
knots for the $x_1$- and $x_2$- dimensions, respectively. It is worth noting that the function $f(x_1, 0)$ is the function considered by Pensoneault et al (2020) which is bounded between 0 and 1.

**Example 5** (Monotonicity in two-dimensional case with large datasets). In this section, we demonstrate the effectiveness of the proposed methodology for handling multidimensional input spaces while processing large datasets. To illustrate this, we use the monotone non-decreasing function $f(x_1, x_2) = 0.9x_1 + \arctan(12x_2)$. This function is particularly interesting because it is linear with respect to the first variable and exhibits greater variation with respect to the second.

<figure><img src="image_9.png" /><figcaption>**Fig. 7** Monotone function estimation for the two-dimensional case. Top: side-by-side comparison of the two estimators, mAP and MAP, with the true function. Bottom: the three top functions at $x_2 = 0$ (left) and $x_1 = 0$ (right). The gray vertical dashed lines represent the locations of the non-uniform one-dimensional knots.</figcaption></figure>

As for boundedness, the target monotone function $f$ is evaluated on a Latin-hypercube design at 900 locations. The training data $\mathbf{y}$ are generated from (1) using the target function $f$ and a noise standard deviation of 0.01% of the range of $\mathbf{y}$. In the top panels of Figure 7, we present side-by-side comparisons of the two proposed estimators, mAP and MAP, along with the true monotone function. In the bottom panels, we illustrate the three top functions at $x_2 = 0$ (left) and $x_1 = 0$ (right). The gray vertical dashed lines indicate the locations of the non-uniform knots. Notably, the proposed approach effectively captures the linearity with respect to the first variable. According to the notation introduction in Section 4, we get $N_1 = 2$ and $s^{(1)} = s^* = \{0, 1\}$. Furthermore, the algorithm prioritizes finding optimal non-uniform knots in the second dimension, which exhibits greater variation. In that case, we found $N_2 = 5$ optimal non-uniform knots (gray vertical dashed lines in the bottom-right panel), with $s^{(2)} = \{0, 0.053, 0.158, 0.421, 1\}$.
<figure><img src="image_10.png" /><figcaption>Fig. 8 The run-time in seconds as a function of the number of samples for generating 6,000 MCMC monotonic surfaces, where the first 1,000 are discarded as burn-in. The number of training samples varies from 900 to 100,000.</figcaption></figure>

Figure 8 illustrates the performance of the proposed approach for handling large datasets. The run-time, measured in seconds, for generating 6,000 MCMC monotonic surfaces is presented as a function of the number of training samples. The proposed method is capable of generating these surfaces with 100,000 training samples in under three minutes, including parameter selection. It is worth noting that the computational run-time increases linearly with the number of observations. Other approaches are not represented because they become prohibitive for a large number of observations $n = 100,000$.

**Example 6** (Monotonicity in 3D with respect to a subset of $\mathcal{X}$). We now investigate the performance of the proposed methodology when monotonicity constraints are imposed on a subset of the input space $\mathcal{X} = [0, 1]^3$ (see Proposition 1). To do this, we consider the following function:

$$f(x_1, x_2, x_3) = 3\{1 + \exp(-10x_1 + 2.1)\} + \cos(5x_2) + x_3.$$

This function is monotone (non-decreasing) with respect to $x_1$ and $x_3$. Additionally, it exhibits more variation with respect to $x_1$ and $x_2$. We expect the proposed approach to firstly activate the input variable $x_2$ and then $x_1$. Additionally, as $f$ is linear with respect to $x_3$, we expect the proposed algorithm to impose only $N_3 = 2$ for the third input variable.

In Figure 9, we illustrate the performance of the proposed approach in the three-dimensional case, where monotonicity constraints are imposed only on $x_1$ and $x_3$. In the top panel, we show side-by-side the true function and the two developed estimators, mAP and MAP, at $\mathbf{x} = (x_1, x_2, 0)$, for any $(x_1, x_2) \in [0, 1]^2$. In the bottom panel, we illustrate these two estimators alongside the true function at $(x_1, 0, 0)$ (left), $(0, x_2, 0)$ (middle), and $(0, 0, x_3)$ (right). The gray vertical dashed lines represent the locations of the one-dimensional knots for each input variable. The proposed methodology effectively captures the monotonicity constraints with respect to $x_1$ and $x_3$, as well as the linear behavior of $f$ with respect to $x_3$. Once again, the proposed algorithm
<figure><img src="image_11.png" /><figcaption>Fig. 9 Monotone function estimation for the three-dimensional case with respect to $x_1$ and $x_3$. Top: side-by-side comparison of the two estimators, mAP and MAP, with the true function at $\boldsymbol{x} = (x_1, x_2, 0)$. Bottom: the proposed estimators at $x_1$ (left), $x_2$ (middle) and $x_3$ (right). The gray vertical dashed lines represent the locations of the non-uniform one-dimensional knots.</figcaption></figure>

places one-dimensional knots at the locations where the function changes behavior. According to Proposition 1, this results in $m = 154$ linear inequality constraints on the Gaussian vector $\boldsymbol{\xi}_S$. **Example 7** (Toy example in higher-dimensions). In this numerical example, the effectiveness of the proposed methodology in high-dimensional input spaces is investigated. To do this, we consider the following 7D target function

$$f(\boldsymbol{x}) = 3/[1 + \exp(-10x_1 + 2.1)] + \arctan(5x_2) + \sum_{i=3}^{7} (i - 1)x_i,$$

where $\boldsymbol{x} = (x_1, \dots, x_7) \in [0, 1]^7$. This function is monotone (non-decreasing) with respect to all its input variables $x_i$, $i = 1, \dots, 7$. It is worth noting that the methodology developed in this paper is capable of incorporating shape constraints with respect to a subset of $\mathcal{X}$ (Example 6). While some computations could potentially be simplified, considering $f$ as an additive function, we choose not to take advantage of this in order to test our approach on a fully 7D example. The target function $f$ is evaluated on a Latin-hypercube design at 5000 locations. The training data $\mathbf{y}$ are generated from (1) using the target function $f$ and a noise standard deviation $\sigma = 0.1$. The number of knots should be chosen to balance high-resolution approximation with minimal computational cost. Based on the algorithm developed in this paper, we found the optimal number of knots per dimension (4, 5, 2, \dots, 2). In a previous work (Maatouk et al, 2025b), we have fixed the number of knots as follows: “we set the number of knots per dimension to two, except for the first two dimensions (i.e., $N_\zeta = 2$, for $\zeta = 3, \dots, 7$), owing to the linearity of the target function concerning the corresponding variables. For the first dimension, we choose $N_1 = 5$, and for the second dimension, we choose $N_2 = 4$, since the function $f$ exhibits more
variation along these dimensions compared to the others". In Maatouk et al (2025b), we have used a uniform discretization (i.e., equally spaced knots per dimension) and thus to a 'rectangle' design knots. However, when the methodology developed in the present paper is considered, we get the following result (Table 1), which are comparable to those obtain in Maatouk et al (2025b).

<table><caption>**Table 1** Comparison of regression error measures for the mAP and MAP estimators using the proposed low-rank approach (left) and Maatouk et al (2025b) (right).</caption><thead><tr><th rowspan="2"></th><th colspan="2">Low-rank</th><th colspan="2">Maatouk et al (2025b)</th></tr><tr><th>mAP</th><th>MAP</th><th>mAP</th><th>MAP</th></tr></thead><tbody><tr><td>MAE</td><td>4.81 × 10<sup>-2</sup></td><td><strong>4.51 × 10<sup>-2</sup></strong></td><td>6.04 × 10<sup>-2</sup></td><td>5.83 × 10<sup>-2</sup></td></tr><tr><td>MSE</td><td>0.36 × 10<sup>-2</sup></td><td><strong>0.31 × 10<sup>-2</sup></strong></td><td>0.65 × 10<sup>-2</sup></td><td>0.57 × 10<sup>-2</sup></td></tr><tr><td>RMSE</td><td>5.98 × 10<sup>-2</sup></td><td><strong>5.59 × 10<sup>-2</sup></strong></td><td>8.09 × 10<sup>-2</sup></td><td>7.56 × 10<sup>-2</sup></td></tr><tr><td>SMSE</td><td>3.99 × 10<sup>-4</sup></td><td><strong>3.47 × 10<sup>-4</sup></strong></td><td>7.91 × 10<sup>-4</sup></td><td>6.90 × 10<sup>-4</sup></td></tr><tr><td>Q<sup>2</sup></td><td>99.96 × 10<sup>-2</sup></td><td><strong>99.97 × 10<sup>-2</sup></strong></td><td>99.92 × 10<sup>-2</sup></td><td>99.93 × 10<sup>-2</sup></td></tr><tr><td>PVA</td><td>1.31</td><td>1.45</td><td><strong>0.90</strong></td><td>1.04</td></tr><tr><td>WAIC</td><td><strong>-7.35 × 10<sup>3</sup></strong></td><td><strong>-7.35 × 10<sup>3</sup></strong></td><td>-6.38 × 10<sup>3</sup></td><td>-6.38 × 10<sup>3</sup></td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Fig. 10** From left to right: the target function, the MAP from Definition 2, and the mAP estimator based on 6,000 MCMC iterations (Definition 1). All these functions are evaluated at $(x_1, x_2, 0, 0, 0, 0)$ (top) and at $(x_1, 0, x_3, 0, \dots, 0)$ (bottom).</figcaption></figure>

Figure 10 displays side-by-side the proposed MAP and mAP estimators together with the true monotone function evaluated at $(x_1, x_2, 0, 0, 0, 0)$. The mAP estimator is obtained through 6,000 MCMC iterations of the algorithm developed in Subsection 4.2, where the first 1,000 iterations are discarded as burn-in. It is worth noting that the proposed approach (8) is capable of capturing the monotonicity dynamics of the target function across the entire domain.
# 6 Real-world data studies

## 6.1 Application to Duncan's Occupational Prestige Data

In this section, the methodology developed in the present paper is applied to Duncan's Occupational Prestige Data (Duncan, 1961). The data for this study are available in the R-package `carData`. The dataset consists of the Prestige scores and other characteristics (Income and Education) of 45 U.S. occupations in 1950. Income represents the percentage of occupational incumbents in the 1950 U.S. Census who earned $3,500 or more per year (approximately $36,000 in 2017 U.S. dollars). Education represents the percentage of occupational incumbents in 1950 who were high school graduates. Finally, Prestige is measured as the percentage of respondents in a social survey who rated the occupation as “good” or better in terms of Prestige. One might reasonably expect that Prestige increases with both Education and Income, and exploratory data analysis confirms this.

<figure><img src="image_13.png" /><figcaption>**Fig. 11** Top: The Prestige as a function of Income (left) and Education (right). Bottom: The mAP estimate (left) and the MAP estimate (right) in the presence of the data (red stars).</figcaption></figure>

In Figure 11, we show the real-world Duncan's Occupational Prestige Data as a function of the two characteristics (Income and Education) and we plot the two competitors estimates mAP (left) and MAP (right) with respect to the two input variables using the proposed approach developed in the present paper.
After normalizing the data (i.e. transforming the input variables to lie within the unit interval), the proposed algorithm begins by activating the second variable (Education) and identifies $N_2 = 7$ optimal, non-uniform, one-dimensional knots. It then proceeds to find $N_1 = 6$ optimal, non-uniform, one-dimensional knots with respect to the first variable (Income). We visually observe that the MAP estimator is more flexible and provides a better fit to the data than the mAP estimator. To quantitatively compare the performance of these two estimators (mAP and MAP), we compute the energy measure criterion:

$$E_{\text{mAP}} = \frac{\sum_{i=1}^{45} [y_i - \bar{Y}_{\text{mAP}}^{S}(\boldsymbol{x}_i)]^2}{\sum_{i=1}^{45} y_i^2} \quad \text{and} \quad E_{\text{MAP}} = \frac{\sum_{i=1}^{45} [y_i - \hat{Y}_{\text{MAP}}^{S}(\boldsymbol{x}_i)]^2}{\sum_{i=1}^{45} y_i^2}.$$

Once again, the MAP estimator slightly outperforms the mAP estimator in terms of prediction accuracy, achieving an energy measure of 0.048 compared to 0.05. It is important to note that the proposed approach can incorporate additional boundedness constraints, in addition to monotonicity, to ensure that the Prestige remains within the interval [0, 100].

## 6.2 Application to coastal flooding (5 dimensional case)

In this section, the methodology developed in the present paper is applied to a 5D coastal flooding case study, analyzed by Azzimonti et al (2019). The data for this study are available in the R-package ProfExtrema (Azzimonti, 2018). Coastal flooding caused by overflow is examined in the Boucholeurs area, located along the French Atlantic coast (Rohmer et al, 2018; Sire et al, 2023). This region, situated near La Rochelle (France) experienced severe flooding during the 2010 Xynthia storm event, resulting in significant loss of life and extensive material damage.

The dataset contains 200 observations of the flooded area ($A_{\text{flood}}$, in square meters (m$^2$)) driven by five offshore forcing conditions (inputs), $x = (T, S, t_0, t_-, t_+)$ at the Boucholeurs area:

* The tide is represented as a sinusoidal signal, parameterized by its high-tide level, $T \in [0.95 \text{ m}, 3.70 \text{ m}]$.

* The surge signal is modeled as a triangular waveform defined by four parameters:

- the peak amplitude $S \in [0.65 \text{ m}, 2.50 \text{ m}]$,

- the phase difference ($t_0$, hours) between the surge peak and the high tide,

- the time duration of the raising part ($t_-$, hours) of the (triangular) surge signal,

- the time duration of the falling part ($t_+$, hours) of the (triangular) surge signal.

It is worth noting that $A_{\text{flood}}$ is known to increase with increasing $T$ and $S$. While the contributions of $T$, $S$, $t_-$, and $t_+$ to $A_{\text{flood}}$ are nearly linear, the variable $t_0$ introduces greater
variability in $A_{flood}$. Based on the methodology developed in this paper, we expect the proposed algorithm to begin discretizing along the variable $t_0$ and to propose more one-dimensional knots in this dimension compared to the others.

<table><caption>**Table 2** Energy measure criterion for the coastal flooding application. Different covariance kernels have been employed.</caption><thead><tr><th rowspan="2">Kernels</th><th>Number of knots</th><th colspan="2">Energy measure criterion</th></tr><tr><th>per dimension</th><th>MAP</th><th>mAP</th></tr></thead><tbody><tr><td>Squared exponential</td><td>(2, 2, 4, 4, 2)</td><td>$E_{MAP} = 4.51 \times 10^{-3}$</td><td>$E_{MAP} = 5.21 \times 10^{-3}$</td></tr><tr><td>Matérn with $\nu = 5/2$</td><td>(2, 2, 5, 4, 2)</td><td>$E_{MAP} = 4.68 \times 10^{-3}$</td><td>$E_{MAP} = 5.44 \times 10^{-3}$</td></tr><tr><td>Matérn with $\nu = 3/2$</td><td>(2, 2, 5, 4, 2)</td><td>$E_{MAP} = 4.76 \times 10^{-3}$</td><td>$E_{MAP} = 5.13 \times 10^{-3}$</td></tr><tr><td>Matérn with $\nu = 1/2$</td><td>(7, 2, 12, 2, 2)</td><td>$E_{MAP} = 3.93 \times 10^{-3}$</td><td>$E_{MAP} = 4.72 \times 10^{-3}$</td></tr></tbody></table>

Table 2 shows the energy measure criterion (Bay et al, 2017) for the two proposed estimators $\bar{Y}_{MAP}^S$ and $\hat{Y}_{MAP}^S$ (mAP and MAP, respectively) together with the number of knots per dimension for different type of covariance functions. Since the target function in this real application $f := \log_{10}(A_{flood})$ is unknown, the energy measure is computed over the 200 available data:

$$E_{MAP} = \frac{\sum_{i=1}^{200} [y_i - \hat{Y}_{MAP}^S(\mathbf{x}_i)]^2}{\sum_{i=1}^{200} y_i^2} \quad \text{and} \quad E_{mAP} = \frac{\sum_{i=1}^{200} [y_i - \bar{Y}_{mAP}^S(\mathbf{x}_i)]^2}{\sum_{i=1}^{200} y_i^2}.$$

As expected, the proposed algorithm adds more one-dimensional knots along the variable $t_0$ and effectively captures the linearity of the other variables. For instance, when using the Squared Exponential (SE) covariance kernel, we obtain 128 knots and basis functions. This result opens the discussion to the following question: does lower regularity of kernels require more knots? This is a particularly interesting result for quantifying uncertainty, as the proposed MCMC algorithm performs better in this context. Notably, the MAP estimator consistently outperforms the mAP estimator in terms of prediction accuracy across all scenarios. This result is coherent with the Bayesian analysis obtained in Maatouk et al (2025a). This performance is attributed to the efficient approach developed in this paper, as well as the parameter selection process, which takes into account the shape constraints. Based on the results in Table 2, the Matérn covariance function with smoothness parameter $\nu = 1/2$ provides more accuracy results compared to the other smoothness covariance kernels. It is worth mentioning that the MAP is an efficient estimator, computable in a couple of seconds since it is based on a QP, whereas the mAP estimator relies on the average of the proposed MCMC algorithm.
# 7 Conclusion

In this paper, multidimensional Bayesian shape-restricted function estimation with large datasets is investigated. Specifically, finite-dimensional Gaussian processes (GPs) are considered, which satisfy a wide range of linear inequality constraints throughout the entire domain. We develop a new methodology that combines a low-rank approach with ordering of the active input variables using the maximum a posteriori estimator. More precisely, we identify an ordering of influential input variables and iteratively compute optimal non-uniform one-dimensional knots. Additionally, we propose an efficient Markov chain Monte Carlo (MCMC) approach based on Elliptical Slice Sampling (ESS), and we estimate the model hyperparameters while taking shape constraints into account. The proposed approach significantly reduces the complexity of constrained GP approximation in multidimensional input settings, enabling it to handle large datasets and address complex shape constraints, even in extreme cases. Moreover, it achieves high prediction accuracy with improved coverage intervals. Its efficiency is demonstrated through experiments on both synthetic and real-world datasets.

# Appendix A Technical proofs

This appendix contains the technical proofs of the lemmas and the corollary presented in the main text.

*Proof of Lemma 1.* For any $x \in \mathcal{X}$, we have

$$ 
\begin{aligned}
|\tilde{f}_N(x) - f(x)(x)| &= \left| \sum_{j=1}^{N} (f(t_j) - f(x))\phi_j(x) \right| \\
&\le \sum_{j=1}^{N} |f(t_j) - f(x)|\phi_j(x) \le \sup_{|x-x'| \le \delta_N} |f(x) - f(x')|.
\end{aligned}
$$

The last inequality holds for the facts $\phi_j(x) = 0$ for $|t_j - x| > \delta_N$ and $\sum_{j=1}^N \phi_j(x) = 1$, for any $x \in \mathcal{X}$. The uniform continuity of the function $f$ on the compact interval $\mathcal{X}$ concludes the proof. □

*Proof of Lemma 2.* Suppose that $k$ is positive-definite. This implies that for any $\tilde{a} \in \mathbb{R}^N$, we have

$$ \tilde{a}^{\top} \tau^2 K \tilde{a} = \tau^2 \sum_{i=1}^{N} \sum_{j=1}^{N} \tilde{a}_i \tilde{a}_j K_{i,j} \ge 0, $$

where $\tilde{a} = [\tilde{a}_1, \dots, \tilde{a}_N]^\top$ and $K_{i,j}$ as defined in (6). This also implies that the eigenvalues of $K$ are non-negative. Let us denote by $K^N \in \mathbb{R}^{J \times J}$ the covariance matrix such that $K_{i,j}^N =$
$k^N(x_i, x_j)$, for any $[x_1, \dots, x_J]^\top \in \mathbb{R}^J$. For any $\tilde{\mathbf{b}} = [\tilde{b}_1, \dots, \tilde{b}_J]^\top \in \mathbb{R}^J$, we have

$$ \tilde{\mathbf{b}}^{\top} \mathbf{K}^{N} \tilde{\mathbf{b}} = \tau^2 \sum_{i=1}^{J} \sum_{j=1}^{J} \tilde{b}_i \tilde{b}_j \phi(x_i)^{\top} \mathbf{K} \phi(x_j) = \tau^2 \tilde{\mathbf{b}}^{\top} \mathbf{K} \tilde{\mathbf{b}} \geq 0, $$

where $\tilde{\mathbf{b}} = \phi \tilde{\mathbf{b}} \in \mathbb{R}^N$, with $\phi \in \mathbb{R}^{N \times J}$ such that $\phi_{,i} = \phi(x_i) = [\phi_1(x_i), \dots, \phi_N(x_i)]^\top$. $\square$

*Proof of Lemma 3.* For any $x \in \mathcal{X}$, we have

$$ \begin{aligned} |\tilde{\mu}_N(x) - \mu(x)| &= |E[\tilde{Y}_N(x)] - E[Y(x)]| = \left| \sum_{j=1}^N E[Y(t_j)]\phi_j(x) - \mu(x) \right| \\ &= \left| \sum_{j=1}^N (\mu(t_j) - \mu(x))\phi_j(x) \right| \\ &\leq \sum_{j=1}^N |\mu(t_j) - \mu(x)|\phi_j(x) \leq \sup_{|x-x'| \leq \delta_N} |\mu(x) - \mu(x')|. \end{aligned} $$

The last inequality holds for the facts $\phi_j(x) = 0$ for $|t_j - x| > \delta_N$ and $\sum_{j=1}^N \phi_j(x) = 1$, for any $x \in \mathcal{X}$. By the uniform continuity of the function $\mu(\cdot)$ on the compact interval $\mathcal{X}$, we deduce that

$$ \sup_{x \in \mathcal{X}} |\tilde{\mu}_N(x) - \mu(x)| \xrightarrow{N \to +\infty} 0. $$

This concludes the first part of the proof. Now, if $\mu$ is Lipchitz, then there exists some finite constant $M > 0$ such that $|\mu(x) - \mu(y)| \leq M|x - y|$, for any $x, y \in \mathcal{X}$. In that case, we have

$$ |\tilde{\mu}_N(x) - \mu(x)| \leq M \sum_{j=1}^{N} |x - t_j| \phi_j(x) \leq M \delta_N = M/(N-1). \quad (A1) $$

The result in (A1) holds uniformly for any $x \in \mathcal{X}$, therefore we complete the proof. $\square$

*Proof of Lemma 4.* Let $\zeta \in [1, d]$ with $d \geq 2$, and let $N_\zeta \geq 2$. For any $x \in [t_1^\zeta, t_{N_\zeta}^\zeta]$, there exists $j \in [1, N_\zeta - 1]$ such that $x \in [t_j^\zeta, t_{j+1}^\zeta]$. In this case, we have $\phi_{j\zeta}^\zeta(x) = 0$ for all $j_\zeta \notin \{j, j+1\}$, and is nonzero otherwise. As a result,

$$ \sum_{j_\zeta=1}^{N_\zeta} \phi_{j\zeta}^{\zeta}(x) = \phi_{j}^{\zeta}(x) + \phi_{j+1}^{\zeta}(x) = \frac{t_{j+1}^{\zeta} - x}{t_{j+1}^{\zeta} - t_{j}^{\zeta}} + \frac{x - t_{j}^{\zeta}}{t_{j+1}^{\zeta} - t_{j}^{\zeta}} = 1. $$

The second equality in (11) is a consequence of (10), which completes the proof. $\square$
Proof of Lemma 5. The proof is similar to that of Lemma 3. The assumption $\max_{j_1} |t_{j_1}^1 - t_{j_1+1}^1| \le \kappa/(N_1 - 1)$ completes the argument. $\square$

Proof of Corollary 1. We consider the function approximation

$$ \tilde{f}_S(\mathbf{x}) = \sum_{j_1=1}^{N_1} \cdots \sum_{j_d=1}^{N_d} f(t_{j_1}^1, \dots, t_{j_d}^d) [\phi_{j_1}^1(x_1) \times \cdots \times \phi_{j_d}^d(x_d)], $$

for $\mathbf{x} = (x_1, \dots, x_d) \in \mathcal{X}$, where each $\phi_{j_\zeta}^\zeta$ is an asymmetric hat basis function defined on a grid $\{t_{j_\zeta}^\zeta\}_{j_\zeta=1}^{N_\zeta}$ over the unit interval, and the function $f$ is continuous on the domain $\mathcal{X} = [0, 1]^d$. To simplify the notation, we define $\mathbf{j} = (j_1, \dots, j_d)$, $\mathbf{t}_j = (t_{j_1}^1, \dots, t_{j_d}^d)$, and

$$ \phi_j(\mathbf{x}) := \prod_{\zeta=1}^{d} \phi_{j_\zeta}^\zeta(x_\zeta). $$

For any $\mathbf{x} \in \mathcal{X}$, by the triangle inequality and (11), we obtain

$$ \begin{aligned} |\tilde{f}_S(\mathbf{x}) - f(\mathbf{x})| &= \left| \sum_j f(t_j) \phi_j(\mathbf{x}) - f(\mathbf{x}) \right| = \left| \sum_j (f(t_j) - f(\mathbf{x})) \phi_j(\mathbf{x}) \right| \\ &\le \sum_j |f(t_j) - f(\mathbf{x})| \phi_j(\mathbf{x}). \end{aligned} $$

Since $\phi_j$ is nonnegative and has compact support near $t_j$, only finitely many terms in the sum are nonzero for each fixed $\mathbf{x}$. Also, because $f$ is uniformly continuous, for every $\varepsilon > 0$, there exists $\delta > 0$ such that

$$ \|\mathbf{x} - t_j\| < \delta \Rightarrow |f(t_j) - f(\mathbf{x})| < \varepsilon. $$

For sufficiently small mesh sizes (i.e., all grid spacings less than $\delta$), all the nonzero terms in the sum satisfy the above inequality, so

$$ |\tilde{f}_S(\mathbf{x}) - f(\mathbf{x})| < \varepsilon. $$

Since this holds for all $\mathbf{x} \in \mathcal{X}$, we conclude that

$$ \sup_{\mathbf{x} \in \mathcal{X}} |\tilde{f}_S(\mathbf{x}) - f(\mathbf{x})| \to 0 \quad \text{as} \quad \max_{\zeta} h_\zeta \to 0, $$

where $h_\zeta = \max_{1 \le j_\zeta < N_\zeta} (t_{j_\zeta}^\zeta - t_{j_\zeta+1}^\zeta)$. Hence, $\tilde{f}_S \to f$ uniformly on $\mathcal{X}$.
From approximation theory, we know that if $f \in H^\alpha(\mathcal{X})$, then there exists a constant $C$ such that the interpolation error in $L^2(\mathcal{X})$ satisfies

$$ \| \tilde{f}_S - f \|_{L^2(\mathcal{X})} \le C h^\alpha \| f \|_{H^\alpha(\mathcal{X})}. $$

For any function $\psi \in C^\infty$ (i.e., infinitely differentiable function), we have by the Cauchy-Schawrz inequality

$$ \left| \int_{\mathcal{X}} (\tilde{f}_S - f)(\mathbf{x}) \psi(\mathbf{x}) d\mathbf{x} \right| \le \|\tilde{f}_S - f\|_{L^2(\mathcal{X})} \times \|\psi\|_{L^2(\mathcal{X})}. $$

Since $\|\tilde{f}_S - f\|_{L^2(\mathcal{X})} \to 0$ at rate $O(h^\alpha)$, and $\psi$ is fixed, then

$$ \int_{\mathcal{X}} \tilde{f}_S(\mathbf{x}) \psi(\mathbf{x}) d\mathbf{x} \to \int_{\mathcal{X}} f(\mathbf{x}) \psi(\mathbf{x}) d\mathbf{x}, $$

so $\tilde{f}_S \to f$ weakly in $L^2(\mathcal{X})$. $\square$

*Proof of Lemma 6.* We have

$$ \begin{aligned} \tilde{\mu} &= E[\xi | \mathbf{X}, \mathbf{y}] = E[\mu] + \mathrm{Cov}(\xi, \mathbf{X}\xi + \epsilon)^{\top} \mathrm{Cov}(\mathbf{X}\xi + \epsilon)^{-1} \mathbf{y} \\ &= \tau^2 (\mathbf{X}\mathbf{K})^{\top} (\tau^2 \mathbf{X}\mathbf{K}\mathbf{X}^{\top} + \sigma^2 \mathrm{I}_n)^{-1} \mathbf{y}. \end{aligned} $$

Additionally, we have

$$ \begin{aligned} \tilde{\Sigma} &= \mathrm{Cov}(\xi | \mathbf{X}, \mathbf{y}) = \mathrm{Cov}(\xi) - \mathrm{Cov}(\xi, \mathbf{X}\xi + \epsilon)^{\top} \mathrm{Cov}(\mathbf{X}\xi + \epsilon)^{-1} \mathrm{Cov}(\xi, \mathbf{X}\xi + \epsilon) \\ &= \tau^2 \mathbf{K} - \tau^4 (\mathbf{X}\mathbf{K})^{\top} (\tau^2 \mathbf{X}\mathbf{K}\mathbf{X}^{\top} + \sigma^2 \mathrm{I}_n)^{-1} (\mathbf{X}\mathbf{K}), \end{aligned} $$

which completes the proof of the Lemma. $\square$

*Proof of Lemma 7.* The proof is based on the following: suppose that $\xi \in \mathbb{R}$ is bounded between $\underline{b} \in \mathbb{R}$ and $\bar{b} \in \mathbb{R}$, i.e., $\underline{b} \le \xi \le \bar{b}$, with $\underline{b} < \bar{b}$ representing lower and upper bounds, respectively. Then,

$$ \begin{aligned} \mathbf{1}_{\{\underline{b} \le \xi \le \bar{b}\}} &\approx \frac{1}{1 + \exp[-\eta(\xi - \underline{b})]} \times \frac{1}{1 + \exp[-\eta(\bar{b} - \xi)]} \\ &= \frac{1}{1 + \exp[-\eta(\xi - \underline{b})] + \exp[-\eta(\bar{b} - \xi)] + \exp[-\eta(\bar{b} - \underline{b})]}, \end{aligned} $$

for a given approximation parameter $\eta > 0$ and for any $\xi \in \mathbb{R}$. $\square$
References

Abrahamsen P, Benth F (2001) Kriging with inequality constraints. Mathematical Geology
33:719–744

Agrell C (2019) Gaussian processes with linear operator inequality constraints. Journal of Machine Learning Research 20(135):1–36

Aronszajn N (1950) Theory of reproducing kernels. Transactions of the American Mathematical Society 68

Azzimonti D (2018) Compute and visualize profile extrema functions. https://cran.r-project.org/web/packages/profExtrema/index.html

Azzimonti D, Ginsbourger D, Rohmer J, et al (2019) Profile extrema for visualizing and quantifying uncertainties on excursion regions: application to coastal flooding. Technometrics

Bachoc F (2013) Cross validation and maximum likelihood estimations of hyper-parameters of Gaussian processes with model misspecification. Computational Statistics & Data Analysis 66:55-69

Bachoc F, López-Lopera A, Roustant O (2022) Sequential construction and dimension reduction of Gaussian processes under inequality constraints. SIAM Journal on Mathematics of Data Science 4(2):772–800

Banerjee A, Dunson D, Tokdar S (2013) Efficient Gaussian process regression for large datasets. Biometrika 100(1):75–89

Bay X, Grammont L, Maatouk H (2016) Generalization of the Kimeldorf-Wahba correspondence for constrained interpolation. Electronic Journal of Statistics 10(1):1580–1595

Bay X, Grammont L, Maatouk H (2017) A new method for interpolating in a convex subset of a Hilbert space. Computational Optimization and Applications 68:95–120

Benavoli A, Azzimonti D (2024) Linearly constrained Gaussian processes are SkewGPs: application to monotonic preference learning and desirability. In: The 40th Conference on Uncertainty in Artificial Intelligence

Berlinet A, Thomas-Agnan C (2011) Reproducing kernel Hilbert spaces in probability and statistics. Springer Science & Business Media
Binois M, Gramacy R (2021) hetgp: Heteroskedastic Gaussian process modeling and sequential design in R. Journal of Statistical Software 98:1–44

Binois M, Gramacy R, Ludkovski M (2018) Practical heteroscedastic Gaussian process modeling for large simulation experiments. Journal of Computational and Graphical Statistics 27(4):808–821

Blanchet-Scalliet C, Demory B, Gonon T, et al (2023) Gaussian process regression on nested spaces. SIAM/ASA Journal on Uncertainty Quantification 11(2):426–451

Botev Z (2017) The normal law under linear restrictions: simulation and estimation via minimax tilting. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 79(1):125–148

Boyd S, Vandenberghe L (2004) Convex optimization. Cambridge University Press

Cai B, Dunson D (2007) Bayesian multivariate isotonic regression splines: applications to carcinogenicity studies. Journal of the American Statistical Association 102(480):1158–1171

Cao J, Katzfuss M (2023) Linear-cost Vecchia approximation of multivariate normal probabilities. arXiv preprint arXiv:231109426

Chataigner M, Cousin A, Crépey S, et al (2021) Beyond surrogate modeling: Learning the local volatility via shape constraints. SIAM Journal on Financial Mathematics 12(3):SC58–SC69

Chen Y, Hosseini B, Owhadi H, et al (2025) Gaussian measures conditioned on nonlinear observations: consistency, MAP estimators, and simulation. Statistics and Computing 35(1):1–23

Chevalier C, Bect J, Ginsbourger D, et al (2014) Fast parallel kriging-based stepwise uncertainty reduction with application to the identification of an excursion set. Technometrics 56(4):455–465

Chipman H, George E, McCulloch R, et al (2022) mBART: Multidimensional Monotone BART. Bayesian Analysis 17(2):515–544

Cotter S, Roberts G, Stuart A, et al (2013) MCMC methods for functions: modifying old algorithms to make them faster. Statistical Science 28(3):424–446
Cousin A, Maatouk H, Rullière D (2016) Kriging of financial term-structures. European Journal of Operational Research 255(2):631-648

Curtis S, Ghosh S (2011) A variable selection approach to monotonic regression with Bernstein polynomials. Journal of Applied Statistics 38(5):961-976

Duncan O (1961) A socioeconomic index for all occupations. In Reiss, A J, Jr (Ed) Occupations and Social Status

Dupuy D, Helbert C, Franco J (2015) DiceDesign and DiceEval: Two R Packages for Design and Analysis of Computer Experiments. Journal of Statistical Software 65(11):1-38

Eilers P, Marx B (1996) Flexible smoothing with B-splines and penalties. Statistical Science 11(2):89-121

Finocchio G, Schmidt-Hieber J (2023) Posterior contraction for deep Gaussian process priors. Journal of Machine Learning Research 24(66):1-49

Genton M (2001) Classes of kernels for machine learning: a statistics perspective. Journal of Machine Learning Research 2:299-312

Goldfarb D, Idnani A (1983) A numerically stable dual method for solving strictly convex quadratic programs. Mathematical Programming 27(1):1-33

Gramacy R (2020) Surrogates: Gaussian process modeling, design, and optimization for the applied sciences. Chapman and Hall/CRC

Grammont L, Maatouk H, Bay X (2024) Equivalence between constrained optimal smoothing and Bayesian estimation. Journal of Nonparametric Statistics 37(1):60-81

Green D, Lindbloom J, Gelb A (2025) Complex-valued signal recovery using a generalized Bayesian LASSO. SIAM/ASA Journal on Uncertainty Quantification 13(2):831-861

Hensman J, Fusi N, Lawrence ND (2013) Gaussian processes for big data. In: Conference on Uncertainty in Artificial Intelligence

Jidling C, Wahlström N, Wills A, et al (2017) Linearly constrained Gaussian processes. Advances in neural information processing systems 30
Johnson A, Jiang D (2018) Shape constraints in economics and operations research. Statistical Science 33(4):527 – 546

Kim J, Luettgen C, Paynabar K, et al (2023) Physics-based penalization for hyperparameter estimation in Gaussian process regression. Computers & Chemical Engineering 178:108320

Kimeldorf G, Wahba G (1970) A correspondence between Bayesian estimation on stochastic processes and smoothing by splines. The Annals of Mathematical Statistics pp 495–502

Kleijnen J, Van Beers W (2013) Monotonicity-preserving bootstrapped kriging metamodels for expensive simulations. Journal of the Operational Research Society 64(5):708–717

Kochan D, Yang X (2025) Blending physics with data using an efficient Gaussian process regression with soft inequality and monotonicity constraints. Frontiers in Mechanical Engineering 10:1410190

Lenk P, Choi T (2017) Bayesian analysis of shape-restricted functions using Gaussian process priors. Statistica Sinica 27(1):43–69

Lenk P, Lee J, Han D, et al (2024) Hierarchical Bayesian spectral regression with shape constraints for multi-group data. Computational Statistics & Data Analysis 200:108036

Maatouk H, Bay X (2016) A new rejection sampling method for truncated multivariate Gaussian random variables restricted to convex sets. In: Cools R, Nuyens D (eds) Monte Carlo and Quasi-Monte Carlo Methods. Springer International Publishing, Cham, pp 521–530

Maatouk H, Bay X (2017) Gaussian process emulators for computer experiments with inequality constraints. Mathematical Geosciences 49(5):557–582

Maatouk H, Rullière D, Bay X (2023) Sampling large hyperplane-truncated multivariate normal distributions. Computational Statistics 39:1779–1806

Maatouk H, Rullière D, Bay X (2024) Truncated multivariate normal distribution under nonlinear constraints, working paper or preprint

Maatouk H, Rullière D, Bay X (2025a) Bayesian analysis of constrained Gaussian processes. Bayesian Analysis 20(3):973–1002
Maatouk H, Rullière D, Bay X (2025b) Efficient constrained Gaussian process approximation using elliptical slice sampling. Bayesian Analysis pp 1–32. https://doi.org/10.1214/25-BA1555

Maatouk H, Rullière D, Bay X (2025c) Large-scale constrained Gaussian processes for shape-restricted function estimation. Statistics and Computing 35(7)

Meyer M, Hackstadt A, Hoeting J (2011) Bayesian estimation and inference for generalised partial linear models using shape-restricted splines. Journal of Nonparametric Statistics 23(4):867–884

Mongwe W, Mbuvha R, Marwala T (2025) Learning equity volatility surfaces using sparse Gaussian processes. In: Bayesian Machine Learning in Quantitative Finance: Theory and Practical Applications. Springer Nature Switzerland, Cham, pp 61–86

Montesinos-López J, Capella A, Christen J, et al (2023) Bayesian estimation of fault slip distribution for slow slip events based on an efficient hybrid optimal directional gibbs sampler and its application to the guerrero 2006 event. Mathematical Geosciences 55(6):859–886

Murray I, Adams R, MacKay D (2010) Elliptical slice sampling. In: Proceedings of the thirteenth international conference on artificial intelligence and statistics, JMLR Workshop and Conference Proceedings, pp 541–548

Neal R (1999) Regression and classification using Gaussian process priors. Bayesian statistics 6 pp 475–501

Pakman A, Paninski L (2014) Exact Hamiltonian Monte Carlo for truncated multivariate Gaussians. Journal of Computational and Graphical Statistics 23(2):518–542

Pensoneault A, Yang X, Zhu X (2020) Nonnegativity-enforced Gaussian process regression. Theoretical and applied mechanics letters 10(3):182–187

Petit S (2022) Improved Gaussian process modeling: Application to Bayesian optimization. PhD thesis, Université Paris-Saclay

Ray P, Pati D, Bhattacharya A (2020) Efficient Bayesian shape-restricted function estimation with constrained Gaussian process priors. Statistics and Computing 30(4):839–853

Riihimaki J, Vehtari A (2010) Gaussian processes with monotonicity information. In: AISTATS, JMLR Proceedings, vol 9. JMLR.org, pp 645–652
Rohmer J, Idier D, Paris F, et al (2018) Casting light on forcing and breaching scenarios that lead to marine inundation: Combining numerical simulations with a random-forest classification approach. Environmental Modelling & Software 104:64-80

Shively T, Walker S, Damien P (2011) Nonparametric function estimation subject to monotonicity, convexity and other shape constraints. Journal of Econometrics 161(2):166-181

Sire C, Le Riche R, Rullire D, et al (2023) Quantizing Rare Random Maps: Application to Flooding Visualization. Journal of Computational and Graphical Statistics 32(4):1556-1571

Spiller E, Wolpert R, Tierz P, et al (2023) The zero problem: Gaussian process emulators for range-constrained computer models. SIAM/ASA Journal on Uncertainty Quantification 11(2):540-566

Taylor J, Benjamini Y (2016) RestrictedMVN: multivariate normal restricted by affine constraints. R package version 1

Titsias M (2009) Variational learning of inducing variables in sparse Gaussian processes. In: AIAS, PMLR, pp 567-574

Ustyuzhaninov I, Kazlauskaitė I, Ek C, et al (2020) Monotonic Gaussian process flows. In: International Conference on Artificial Intelligence and Statistics, PMLR, pp 3057-3067

Williams C, Rasmussen C (2006) Gaussian processes for machine learning, vol 2. MIT press
Cambridge, MA

Williams N, Osborne C, Seymour I, et al (2023) Application of finite Gaussian process distribution of relaxation times on SOFC electrodes. Electrochemistry Communications p 107458

Wood A, Chan G (1994) Simulation of stationary Gaussian processes in $[0, 1]^d$. Journal of Computational and Graphical Statistics 3(4):409-432

Zhou S, Giulani P, Piekarewicz J, et al (2019) Reexamining the proton-radius problem using constrained Gaussian processes. Physical Review C 99:055202

Zhou S, Ray P, Pati D, et al (2024) A mass-shifting phenomenon of truncated multivariate normal priors. Journal of the American Statistical Association 119(545):582-596