# Sequential construction and dimension reduction of Gaussian processes under constraints

François Bachoc, Andrés F. López-Lopera, Olivier Roustant

► To cite this version:

François Bachoc, Andrés F. López-Lopera, Olivier Roustant. Sequential construction and dimension reduction of Gaussian processes under constraints. SIAM Journal on Mathematics of Data Science, 2022, 4 (2), pp.772-800. 10.1137/21M1407513. hal-03614303

HAL Id: hal-03614303

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

Submitted on 29 Oct 2024

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Sequential construction and dimension reduction of Gaussian
processes under constraints

François Bachoc*, Andrés F. López Lopera† and Olivier Roustant‡

*Inst. de Mathématiques de Toulouse (IMT), UMR5219 CNRS, 31062 Toulouse, France

†Univ. Polytechnique Hauts-de-France, Céramaths, F-59313 Valenciennes, France

‡IMT, UMR5219 CNRS, INSA, 31077 Toulouse cédex 4, France

March 11, 2022

Abstract

Accounting for inequality constraints, such as boundedness, monotonicity or convexity, is chal-
lenging when modeling costly-to-evaluate black box functions. In this regard, finite-dimensional
Gaussian process (GP) regression models bring a valuable solution, as they guarantee that the
inequality constraints are satisfied everywhere. Nevertheless, these models are currently restricted
to small dimensional situations (up to dimension 5). Addressing this issue, we introduce the Max-
Mod algorithm that sequentially inserts one-dimensional knots or adds active variables, thereby
performing at the same time dimension reduction and efficient knot allocation. We prove the con-
vergence of this algorithm. In intermediary steps of the proof, we propose the notion of multi-affine
extension and study its properties. We also prove the convergence of finite-dimensional GPs, when
the knots are not dense in the input space, extending the recent literature. With simulated and real
data, we demonstrate that the MaxMod algorithm remains efficient in higher dimension (at least in
dimension 20), and needs fewer knots than other constrained GP models from the state-of-the-art,
to reach a given approximation error.

# 1 Introduction

Gaussian processes (GPs) are widely used to address diverse applications since they form a flexible prior over functions [45, 51]. They have been successfully applied in research fields such as numerical code approximations [49], global optimization [32], model calibration [34], geostatistics [12, 43] and machine learning [45].

It is known that accounting for inequality constraints (e.g. positivity, monotonicity, convex-
ity) in GPs leads to smaller prediction errors and to more realistic uncertainties [15, 16, 24, 39,
42, 46, 47]. These inequality constraints correspond to available information on functions over
which GP priors are considered. They are encountered in diverse research fields such as social
system analysis [47], computer networking [24], econometrics [13], geostatistics [42], nuclear safety
criticality assessment [39], tree distributions [40], coastal flooding [38], and nuclear physics [52].
Note also that, beyond GPs, regression methods that account for these inequality constraints, also
called shape constraints, are acknowledged as an important need in statistics and machine learning
[8, 20, 21, 27, 28, 29, 36].

Among the existing approaches which enable us to impose inequality constraints to GP mod-
els, we focus on those based on the approximation of GP samples in finite-dimensional spaces of
functions such as piecewise linear functions [3, 13, 39, 42, 52]. Indeed, the main benefit of these
approaches is that they guarantee the inequality constraints to be satisfied everywhere in the input
space. For instance, in the case of boundedness constraints, this means that the realizations from
the posterior distribution of a constrained GP model obtained from these approaches are above or
below the prescribed bounds everywhere in the input space. In contrast, realizations from the pos-
terior distributions obtained from several other GP-based approaches are guaranteed to be above
or below the prescribed bounds only at a limited number of selected input points (see, e.g., [15, 16]
and the work in [47] accounting for monotonicity constraints). Hence, the approaches based on the
approximation of GP samples in finite-dimensional spaces take into account the full information of
the inequality constraints. Furthermore, in practice, domain experts may consider that a statistical
model on a black box function is more trustworthy and physically interpretable if it respects known
inequality constraints everywhere in the input space. For instance, in the application of Section 5.3, a prediction of the flooded area is expected to be non-negative and non-decreasing with respect to the tide and surge inputs. From this point of view, the approximation of GP samples in finite-dimensional spaces is one of the few admissible methods.

Nevertheless, the main drawback of approximating GP samples in finite-dimensional spaces is the scalability to high-dimensional input spaces. Indeed, the finite-dimensional spaces rely on basis functions, each of them being centered at a *D*-dimensional knot, with *D* being the input space dimension. These *D*-dimensional knots need to be obtained from the tensorization of *D* sets of one-dimensional knots, in order to satisfy the constraints everywhere. According to the state-of-the-art, for instance [38, 39, 42], the sets of knots are fixed a priori for each of the *D* inputs. This limits the applicability to small dimension, say, 3 to 5.

In this paper, we overcome this limitation, in situations where the dimension *D* is allowed to be significantly larger (for instance *D* = 20 in Section 5), but where there are many irrelevant input variables, or, in other words, the effective dimension is small. We suggest a sequential procedure for knot insertion and variable selection that is scalable to these situations of higher dimensional input spaces with irrelevant variables. The procedure leverages three important intuitive principles. First, there should be a higher concentration of knots in input regions where the function is varying most. Second, the most influential variables should be allocated the most knots. Third, weakly influential variables should be allocated one-dimensional knots with the least priority.

Let us now describe the sequential procedure. Consider a set of *n* input points and *n* corresponding observations of the function of interest, to be interpolated, as well as given inequality constraints. We start with a coarse finite-dimensional GP model based on few active variables and small sets of one-dimensional knots for them. Then, at each step, we either add a new active variable or insert a one-dimensional knot to a variable that is already active. The new variable or new knot is the one that corresponds to the largest modification, in $L^2$ norm, of the maximum a posteriori (MAP) function, also called the mode, of the constrained GP model. For this reason the suggested sequential procedure is called the MaxMod (maximum modification) algorithm, hereafter called MaxMod. The MAP is the most probable function that interpolates the observations and that satisfies the inequality constraints everywhere, according to the finite-dimensional constrained GP model. As shown in [42], computing the MAP yields a convex optimization problem of moderate complexity. Here, we provide a computationally simple expression of the subsequent $L^2$ norm (in Appendix B), resulting in a computational complexity that is linear in the number of multi-dimensional knots. We allow for free locations of one-dimensional knots by using asymmetric hat basis functions, instead of the symmetric ones investigated in [39, 42]. The sequential procedure also naturally incorporates a penalization for adding new variables, or for inserting one-dimensional knots, that overly increase the total number of multi-dimensional knots.

From the point of view of free knot insertion in spline approximation, MaxMod differs from many existing references [14, 17, 18, 26, 30, 33, 35, 50] (that typically address spline approximation independently of GPs and inequality constraints). Indeed, these references are based on directly evaluating and minimizing the approximation error of a target function and thus rely on multiple evaluations of this function. In contrast, MaxMod is adapted to the situation where evaluations of the target function are scarce, and it simply maximizes the difference between successive spline approximations.

We provide a convergence guarantee for MaxMod. We consider the set of *n* input points and function observations to be fixed and we let the number of iterations go to infinity. This corresponds to increasing the computational budget, as measured by the number of multi-dimensional knots. Then, we show that all the variables are eventually activated and that the set of *D*-dimensional knots becomes dense in the input space. This implies, based on [5, 6], that the MAP function obtained from MaxMod converges to the optimal constrained interpolant function in the reproducing kernel Hilbert space (RKHS, see for instance [10]) of the covariance function of the GP model. Hence, the convergence result states that MaxMod, despite being a sequential procedure, becomes globally efficient as the number of iterations increases. In particular, loosely speaking, the procedure does not fall into an undesirable local pattern, where the inserted knots would cluster and would not eventually cover the whole input space. It is, in general, important to ensure that sequential procedures avoid undesirable local patterns [7, 9].

In order to obtain the convergence result, we extend the results of [5, 6]. These results show that, given a dense sequence of multi-dimensional knots, the MAP of a constrained finite-dimensional GP converges to the above discussed optimal constrained interpolant function in the RKHS of the GP covariance function. The extension tackles the case where the multi-dimensional knots are not dense. Based on the subset *F* of the input space corresponding to their closure, we define a transformation that we call the multiaffine extension, that extends a function defined on *F* to the
entire input space. This extension enables us to define an optimal constrained interpolant function based on a new RKHS restricted to $F$. Then, we show that the MAP converges to this optimal constrained interpolant function, thus extending [5, 6] to any sequence of multi-dimensional knots, not necessarily dense on the entire input space. The construction and properties of the multiaffine extension and this extension of [5, 6] may be of independent interest. Also, this general proof scheme for the convergence of MaxMod may be adapted to the convergence of other algorithms based on hat basis functions.

The benefit of the suggested sequential procedure is shown in a series of numerical experiments, provided in Section 5, with simulated and real-world data. For the latter, data come from a coastal flooding application (see [2, 38]) satisfying both positivity and monotonicity constraints. We test the versatility of MaxMod for efficiently inserting knots or adding active dimensions while reducing the approximation error of the resulting constrained GP. We demonstrate that MaxMod remains tractable and yields a constrained GP model with accurate predictions, even up to the dimension $D = 20$, for which the state-of-the-art procedures either are intractable [39, 42] or do not satisfy the constraints everywhere [16, 47]. Even in smaller dimension, when the procedures of [39, 42] are tractable, MaxMod typically needs fewer knots to achieve a comparable approximation error. We also show the benefit of having fewer knots when subsequently computing confidence intervals from constrained GP models.

This paper is organized as follows. In Section 2, we describe the finite-dimensional GP approach proposed in [39, 42], that we adapt to the case where only a subset of the $D$ variables is active. In Section 3, we introduce MaxMod. In Section 4, we present the multiaffine extension and establish the various convergence results. The numerical experiments are carried out in Section 5. Section 6 concludes the paper. Finally, in the appendix, we provide technical developments, some of the technical conditions, and all the proofs of the paper.

# 2 Finite-dimensional constrained Gaussian processes

For convenience, we have summarized the notations of Sections 2 and 3 in Table 1, located at the end of Section 3.

## 2.1 Basis function decomposition

The principle of finite-dimensional constrained GPs is to consider linear combinations of basis functions which are tensorizations of one-dimensional asymmetric hat basis functions. These one-dimensional basis functions are parametrized by $-\infty < u < v < w < +\infty$, and are written $\phi_{u,v,w}: \mathbb{R} \to \mathbb{R}$, defined by,

$$ \phi_{u,v,w}(t) = \begin{cases} \frac{1}{v-u}(t-u) & \text{for } u \le t \le v \\ \frac{1}{w-v}(w-t) & \text{for } v \le t \le w \\ 0 & \text{for } t \notin [u, w], \end{cases} $$

for $t \in \mathbb{R}$. Clearly, $\phi_{u,v,w}$ is a ‘hat’ function centered at $v$ and with support $[u, w]$.

Still in dimension one, we now explain how a set of basis functions can be defined from the notion of subdivision. A (one-dimensional) subdivision is a set of (one-dimensional) knots $S = \{t_0^{(S)}, \dots, t_{m_S+1}^{(S)}\}$. We let $t_{(0)}^{(S)} \le \dots \le t_{(m_S+1)}^{(S)}$ be the corresponding ordered knots. For any subdivision, we assume that $0 = t_{(1)}^{(S)} < \dots < t_{(m_S)}^{(S)} = 1$ and, to deal with boundary issues, we set $t_{(0)}^{(S)} = -1$ and $t_{(m_S+1)}^{(S)} = 2$. With this convention, the smallest subdivision containing 0 and 1, denoted as $S^0$, is $\{-1, 0, 1, 2\}$.

Now, let $D \in \mathbb{N}$ be the ambient (potentially large) dimension and consider a given set of active variables $\mathcal{J} = (a_1, \dots, a_d) \subseteq \{1, \dots, D\}$, with $1 \le a_1 < \dots < a_d \le D$, of size $d \le D$. We now explain how to construct $d$-dimensional basis functions of the $d$ active variables. We define a $d$-dimensional subdivision (indexed by $\mathcal{J}$) as a vector of one-dimensional subdivisions $S = (S_{a_1}, \dots, S_{a_d})$. For convenience, we may identify $S$ with the set of $d$-dimensional knots $\prod_{j=1}^d S_{a_j} \subseteq \mathbb{R}^d$. We denote by $S_{\mathcal{J}}$ the set of $d$-dimensional subdivisions indexed by $\mathcal{J}$.

For conciseness, we use tensor notation. Thus, the notation for a multi-index is $\ell = (\ell_{a_1}, \dots, \ell_{a_d}) \in \mathbb{N}^d$. For a subdivision $S \in S_{\mathcal{J}}$, the associated set of multi-indices is denoted

$$ \mathcal{L}_S = \{\ell; \ell_i \in \{1, \dots, m_{S_i}\}, i \in \mathcal{J}\}, $$
<figure><img src="image_2.png" /><figcaption>Figure 1: Illustration of the knots and basis functions.</figcaption></figure>

and $A_S = \{\alpha_{\ell} \in \mathbb{R}; \ell \in \mathcal{L}_S\}$ is the set of associated real-valued sequences. For a given multi-index $\ell \in \mathcal{L}_S$, the associated vector of knots is denoted $t_{\ell}^{(S)} = (t_{(\ell a_1)}^{(S_{a_1})}, \dots, t_{(\ell a_d)}^{(S_{a_d})})$. We call this vector a $d$-dimensional knot. Then, to any $\ell \in \mathcal{L}_S$ we associate the $d$-dimensional basis function defined by tensorization, for $t = (t_{a_1}, \dots, t_{a_d}) \in [0, 1]^d$,

$$ \phi_{\ell}^{(S)}(t) = \prod_{i=1}^{d} \phi_{t_{(\ell a_i-1)}^{(S_{a_i})}, t_{(\ell a_i)}^{(S_{a_i})}, t_{(\ell a_i+1)}^{(S_{a_i})}} (t_{a_i}). $$

This basis function is 1 at the knot $(t_{(\ell a_1)}^{(S_{a_1})}, \dots, t_{(\ell a_d)}^{(S_{a_d})})$ and has support equal to the hypercube $\prod_{i=1}^{d} [t_{(\ell a_i-1)}^{(S_{a_i})}, t_{(\ell a_i+1)}^{(S_{a_i})}]$.

In Figure 1, we now provide an illustrative example. Consider a total of $D = 4$ input variables, among which the $d = 2$ variables 1 and 3 are active. We thus have the active set $\mathcal{J} = (a_1, a_2) = (1, 3)$. Consider that for the first active variable, the subdivision is $S_1 = \{-1, 0, 1/3, 1, 2\}$. In $S_1$, the knots $-1$ and $2$ are only present to deal with boundary issues and the knots $t_{(1)}^{(S_1)} = 0$, $t_{(2)}^{(S_1)} = 1/3$ and $t_{(3)}^{(S_1)} = 1$ are displayed as two blue bullets and one blue triangle on the left panel. The basis functions $\phi_{-1,0,1/3}$, $\phi_{0,1/3,1}$ and $\phi_{1/3,1,2}$ centered at these latter three knots are also displayed. Similarly for the second active variable, the subdivision is $S_3 = \{-1, 0, 1/4, 1/2, 1, 2\}$ (this subdivision is written $S_3$ because the second active variable is the variable 3). The four knots $t_{(1)}^{(S_3)} = 0$, $t_{(2)}^{(S_3)} = 1/4$, $t_{(3)}^{(S_3)} = 1/2$ and $t_{(4)}^{(S_3)} = 1$ are displayed at the middle panel, as well as the basis functions $\phi_{-1,0,1/4}$, $\phi_{0,1/4,1/2}$, $\phi_{1/4,1/2,1}$ and $\phi_{1/2,1,2}$ centered at these knots.

Then, consider the two-dimensional space formed by the two active variables with subdivision $S = (S_1, S_3)$. Its associated set of multi-indices $\mathcal{L}_S$ contains the $3 \times 4 = 12$ indices $(1, 1), (1, 2), \dots, (3, 4)$. The corresponding 12 bi-dimensional knots are displayed on the right. The bi-dimensional knot represented as a triangle is $t_{(2,3)}^{(S)} = (1/3, 1/2)$ with multi-index $\ell = (2, 3)$. It corresponds to the triangle one-dimensional knots $t_{(2)}^{(S_1)} = 1/3$ and $t_{(3)}^{(S_3)} = 1/2$. The bi-dimensional basis function $\phi_{\ell}^{(S)}$ on the right panel is a pyramid which top is located at the blue triangle and which vertices are the dashed lines.

Next, let us construct the finite-dimensional space of functions $E_S$ that contains finite-dimensional constrained GPs. We let $E_S$ be the linear space of functions $[0, 1]^d \to \mathbb{R}$, spanned by the tensor basis functions $\phi_{\ell}^{(S)}$ with $\ell \in \mathcal{L}_S$. Equivalently, $E_S$ is the space of multivariate splines of degree 1, constituted of componentwise piecewise linear functions, with knots defined by $S$. For $\alpha \in A_S$, we let $Y_{S,\alpha}$ be the element of $E_S$ with coefficients $\alpha$:

$$ Y_{S,\alpha} = \sum_{\ell \in \mathcal{L}_S} \alpha_{\ell} \phi_{\ell}^{(S)}. \qquad (1) $$

Note that $Y_{S,\alpha}(t_{\ell}^{(S)}) = \alpha_{\ell}$ for all $\ell \in \mathcal{L}_S$. Finally, for a function $f: [0, 1]^d \to \mathbb{R}$, we denote by $\pi_S(f)$ the projection of $f$ onto $E_S$:

$$ \pi_S(f) = \sum_{\ell \in \mathcal{L}_S} f(t_{\ell}^{(S)}) \phi_{\ell}^{(S)}. \qquad (2) $$

This projection is the componentwise piecewise linear function that coincides with $f$ at the knots of $S$.
## 2.2 Function spaces for interpolation and inequality constraints

We let $\mathcal{F}(A, \mathbb{R})$ (resp. $\mathcal{C}(A, \mathbb{R})$) be the set of functions (resp. continuous functions) from a subset $A$ of a finite-dimensional vector space to $\mathbb{R}$. For $d \in \{1, \dots, D\}$, $U = (u_1, \dots, u_n) \in ([0, 1]^d)^n$ and $v^{(n)} = (v_1, \dots, v_n) \in \mathbb{R}^n$, we write

$$ \mathcal{I}_{U,v^{(n)}} = \{f : [0, 1]^d \to \mathbb{R}; f(u_i) = v_i \text{ for } i = 1, \dots, n\}, $$

the set of $d$-dimensional functions which interpolate $v^{(n)}$ at $U$. Typically, a $D$-dimensional function of interest is observed at observation points $\delta_1, \dots, \delta_n \in [0, 1]^D$, with values $v_1, \dots, v_n$ and we let $u_i = (\delta_i)_J$. Throughout the paper, we use the following notation: for $J \subseteq \{1, \dots, D\}$ and $x \in \mathbb{R}^D$, we let $x_J$ be the vector extracted from $x$ by keeping only the components with indices in $J$.

We treat the inequality constraints as a subset $C_D$ of $\mathcal{C}([0, 1]^D, \mathbb{R})$. Note that, even if a constrained GP model has a subset $J$ of active variables that is strictly smaller than $\{1, \dots, D\}$, the functions of this model are also functions of the full $D$ variables, and considered as such, they are required to belong to $C_D$.

Three classical examples of inequality constraints are given below, corresponding respectively to boundedness, (componentwise) monotonicity and componentwise convexity:

$$ C_D = \{f \in \mathcal{C}([0, 1]^D, \mathbb{R}); a \le f(x) \le b \text{ for all } x \in [0, 1]^D\}, \quad (3) $$

$$ C_D = \{f \in \mathcal{C}([0, 1]^D, \mathbb{R}); f(u) \le f(v) \text{ for all } u, v \in [0, 1]^D, u \le v\}, \quad (4) $$

$$ C_D = \{f \in \mathcal{C}([0, 1]^D, \mathbb{R}); \text{for all } i \in \{1, \dots, D\}, \text{ for all } x_{\sim i} \in [0, 1]^{D-1}, \quad (5) $$

the function $u_i \mapsto f(u_i, x_{\sim i})$ is convex\}.

For monotonicity, the notation used above, $u \le v$, means $u_1 \le v_1, \dots, u_D \le v_D$. For componentwise convexity, for $t \in [0, 1]^D$, $i \in \{1, \dots, D\}$ and $u_i \in [0, 1]$, we denote $t_{\sim i}$ the vector obtained from $t$ by removing the coordinate $i$, and $(u_i, t_{\sim i})$ the vector obtained from $t$ by replacing the $i$th coordinate by $u_i$.

The inequality constraint set $C_D$ of functions in $\mathcal{C}([0, 1]^D, \mathbb{R})$ naturally yields a corresponding inequality constraint set for functions in $\mathcal{C}([0, 1]^{|J|}, \mathbb{R})$, where for a finite set $\Theta$, we write $|\Theta|$ for its cardinality. The inequality constraint set in $\mathcal{C}([0, 1]^{|J|}, \mathbb{R})$ corresponding to $C_D$ is

$$ C_J = \{g : [0, 1]^{|J|} \to \mathbb{R}; (x \in [0, 1]^D \mapsto g(x_J)) \in C_D\}. $$

The set $C_J$ is the set of functions of the variables in $J$, that yield functions in $C_D$ when extended with $D - |J|$ inactive variables.

We finally introduce the next basic condition on the inequality constraint set.

**Condition 1** (constraint set topology). The set $C_D$ is convex and closed with the topology of uniform convergence.

Note that Condition 1 holds when $C_D$ is given by one of (3), (4) or (5).

## 2.3 Finite-dimensional Gaussian processes under interpolation and inequality constraints

### 2.3.1 Finite-dimensional Gaussian processes

We consider a continuous GP $\xi_D$ indexed by $[0, 1]^D$, with mean zero and continuous covariance function $k_D$ on $[0, 1]^D \times [0, 1]^D$. The covariance function $k_D$ can be restricted to an active set of variables $J \subseteq \{1, \dots, D\}$ as follows. We let $-J = \{1, \dots, D\} \setminus J$. We let $k_J$ be the covariance function on $[0, 1]^{|J|} \times [0, 1]^{|J|}$ defined by, for $u, v \in [0, 1]^{|J|}$, $k_J(u, v) = k_D(\tilde{u}, \tilde{v})$, where $\tilde{u}_J = u$, $\tilde{v}_J = v$ and $\tilde{u}_{-J} = 0$, $\tilde{v}_{-J} = 0$ and $\tilde{u}_{-J} = 0$. Note that the choice of the value 0 for the inactive variables, i.e. that are not in $J$, is arbitrary.

We consider the continuous GP $\xi_J$ on $[0, 1]^{|J|}$ defined by, for $u \in [0, 1]^{|J|}$, $\xi_J(u) = \xi_D(\tilde{u})$, where $\tilde{u} \in [0, 1]^D$ is defined by $\tilde{u}_J = u$ and $\tilde{u}_{-J} = 0$. Then $\xi_J$ has mean zero and covariance function $k_J$. Remark that $\xi_J$ is obtained by freezing to 0 the inputs of $\xi_D$ that are not in $J$, which enables to model functions of the variables in $J$. Given a multi-dimensional subdivision, $S \in S_J$, we consider the finite-dimensional GP

$$ \pi_S(\xi_J) = \sum_{\ell \in \mathcal{L}_S} \xi_J(t_{\ell}^{(S)}) \phi_{\ell}^{(S)}. \quad (6) $$

This finite-dimensional GP only depends on the vector of values at the knots, $(\xi_J(t_{\ell}^{(S)}))_{\ell \in \mathcal{L}_S}$.
Let us consider the covariance matrix of this vector, that we write $k_J(S, S)$. Writing $J = (a_1, \dots, a_d)$ with $1 \le a_1 < \dots < a_d \le D$, $k_J(S, S)$ is the matrix of size $m_{S_{a_1}} \times \dots \times m_{S_{a_d}}$, that we write in a multi-index way as, for $\ell, \ell' \in \mathcal{L}_S$, $k_J(S, S)_{\ell, \ell'} = k_J(t_{\ell}^{(S)}, t_{\ell'}^{(S)})$.

Remark that with this multi-index writing, matrix products of the form $k_J(S, S)k_J(S, S)$, matrix inverses of the form $k_J(S, S)^{-1}$ and matrix vector products of the form $k_J(S, S)\alpha$ for $\alpha \in A_S$ can be defined by a straightforward extension of the corresponding operations for standard (single indexed) matrices and vectors. Let us explain the matrix vector product case. Since $\mathcal{L}_S$ corresponds to the set $\prod_{i=1}^d \{1, \dots, m_{S_{a_i}}\}$ and $\alpha = (\alpha_\ell)_{\ell \in \mathcal{L}_S}$, the values of $\alpha$ can be indexed (arbitrarily) as $\alpha_1, \dots, \alpha_N$, with $N = m_{S_{a_1}} \times \dots \times m_{S_{a_d}}$. With the same indexing, the values in $k_J(S, S)$ can be written as $(k_{i,j})_{i,j=1,\dots,N}$. The standard matrix vector product between $(k_{i,j})_{i,j=1,\dots,N}$ and $(\alpha_1, \dots, \alpha_N)$ provides the vector $(\beta_1, \dots, \beta_N)$. Using then the reverse indexing, we obtain the values $(\beta_\ell)_{\ell \in \mathcal{L}_S}$ which are exactly the values of the (multi-indexed) matrix vector product $k_J(S, S)\alpha$.

We assume that $k_J(S, S)$ is invertible for all $J \subseteq \{1, \dots, D\}$ and $S \in \mathcal{S}_J$. This holds in particular when the matrix $[k_D(\delta_i, \delta_j)]_{i,j=1,\dots,q}$ is invertible for any $q \in \mathbb{N}$ and $\delta_1, \dots, \delta_q \in [0, 1]^D$, two-by-two distinct. This is verified by most common covariance functions, for instance the squared exponential ones and those from the Matérn class [45, 51].

### 2.3.2 Obtaining a finite number of linear inequality constraints for finite-dimensional Gaussian processes

The main benefit of the finite-dimensional GP $\pi_S(\xi_J)$ is that, for many classical inequality sets $\mathcal{C}_J$, obtained from $J = (a_1, \dots, a_d) \subseteq \{1, \dots, D\}$ and $\mathcal{C}_D$, there exists an explicit (multi-indexed) matrix $M(\mathcal{C}_J) = (M(\mathcal{C}_J)_{b,\ell})_{b=1,\dots,B}$ and an explicit vector $v(\mathcal{C}_J) = (v(\mathcal{C}_J)_b)_{b=1,\dots,B}$ such that

$$ \pi_S(\xi_J) \in \mathcal{C}_J \iff M(\mathcal{C}_J)(\xi_J(t_{\ell}^{(S)}))_{\ell \in \mathcal{L}_S} \le v(\mathcal{C}_J), \quad (7) $$

where again the definition of multi-indexed matrix-vector products is straightforward.

Equation (7) provides $B$ linear inequality constraints on the vector of values of $\xi_J$ at the knots $(t_{\ell}^{(S)})_{\ell \in \mathcal{L}_S}$. Hence, the constraint $\pi_S(\xi_J) \in \mathcal{C}_J$, that is a priori infinite-dimensional and intractable, boils down to simple linear inequality constraints.

Consider for illustration the one-dimensional case, i.e. $d=1$, and let $S = S_{a_1}$. Then, in the case of boundedness constraints, it is shown in [42] that $\pi_S(\xi_J) \in \mathcal{C}_J$ if and only if

$$ \xi_J(t_{\ell}^{(S)}) \in [a, b], \quad \ell = 1, \dots, m_S. $$

Hence, for boundedness constraints, (7) holds when $M(\mathcal{C}_J)$ is the (single-indexed) $2m_S \times m_S$ matrix $(-I_{m_S}, I_{m_S})^\top$ and $v(\mathcal{C}_J)$ is $(-a, \dots, -a, b, \dots, b)^\top$ of size $2m_S \times 1$.

In the case of monotonicity constraints, it is shown in [42] that $\pi_S(\xi_J) \in \mathcal{C}_J$ if and only if

$$ \xi_J(t_{\ell}^{(S)}) \ge \xi_J(t_{\ell-1}^{(S)}), \quad \ell = 2, \dots, m_S. $$

Hence under monotonicity constraints, (7) holds when $v(\mathcal{C}_J)$ is the zero vector and $M(\mathcal{C}_J)$ is the (single-indexed) $(m_S - 1) \times m_S$ banded matrix

$$ \begin{pmatrix} 1 & -1 & 0 & \cdots & 0 \\ 0 & 1 & -1 & 0 & \cdots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & \cdots & 0 & 1 & -1 \end{pmatrix}. $$

In dimension $|J| = d > 1$, for boundedness and monotonicity constraints, the principle is the same but the notations become more cumbersome. In Appendix A, we provide the expressions of $M(\mathcal{C}_J)$ and $v(\mathcal{C}_J)$ for which (7) holds, for boundedness, monotonicity and componentwise convexity, in any dimension. These expressions follow from [39, 42], except for componentwise convexity in dimension larger than one, when they are not available in earlier references, to the best of our knowledge.
### 2.3.3 The maximum a posteriori function

Consider now some input locations $U = (u_1, \dots, u_n) \in ([0, 1]^d)^n$ and some corresponding observations $v^{(n)} = (v_1, \dots, v_n) \in \mathbb{R}^n$. We are interested in the conditional distribution of $\pi_S(\xi_J)$ given $\pi_S(\xi_J) \in \mathcal{I}_{U,v^{(n)}}$ (interpolation constraints) and given $\pi_S(\xi_J) \in \mathcal{C}_J$ (inequality constraints).

Let us define the following central characteristics of this conditional distribution:

$$ \hat{\alpha}_{J,S,U,v^{(n)}} = \underset{\alpha \in A_S}{\operatorname{argmin}} \alpha^\top k_J(S, S)^{-1}\alpha. \quad (8) $$
$$ \text{s.t.} \quad Y_{S,\alpha} \in \mathcal{I}_{U,v^{(n)}} \cap \mathcal{C}_J $$

The mathematical definitions underlying (8) are recalled in Table 1. In words, the quantity in (8) is the mode of the density of the Gaussian vector extracted from the finite-dimensional GP at the knots, conditionally on the equality and inequality constraints. The linear combination of the basis functions based on these knot values leads to the MAP function

$$ Y_{S,\hat{\alpha}_{J,S,U,v^{(n)}}}, \quad (9) $$

which is used by [5, 6, 39, 42] and is also called the mode of a finite-dimensional GP given interpolation and inequality constraints.

Note that the interpolation constraint $Y_{S,\alpha} \in \mathcal{I}_{U,v^{(n)}}$ can be expressed as a set of explicit linear equations for $\alpha$, see for instance [39, 42]. When the equivalence (7) holds, the constraint $Y_{S,\alpha} \in \mathcal{C}_J$ is expressed as a set of linear inequality constraints. In this case, the optimization problem (8) is a quadratic optimization problem with linear inequality constraints and efficient optimization procedures are available [11, 25], see also [39, 42].

Furthermore, under Condition 1 and when $u_1, \dots, u_n$ are two-by-two distinct, $\hat{\alpha}_{J,S,U,v^{(n)}}$ in (8) is well-defined when the minimization set is non-empty, because this set is closed and convex, and the function to be minimized is continuous and strictly convex, and goes to infinity as $||\alpha||$ goes to infinity.

Finally, beside the MAP function $Y_{S,\hat{\alpha}_{J,S,U,v^{(n)}}}$, obtaining conditional realizations of $\pi_S(\xi_J)$ given $\pi_S(\xi_J) \in \mathcal{I}_{U,v^{(n)}}$ and $\pi_S(\xi_J) \in \mathcal{C}_J$ is also of high interest, for uncertainty quantification. When (7) holds, these conditional realizations can be approximately sampled by Monte Carlo and Markov chain Monte Carlo (MCMC) procedures [39, 42] (see also Section 5).

## 3 The MaxMod algorithm

### 3.1 Initialization

We consider the set of $n$ observation points and observed values to be fixed in Section 3. We let $x_D^{(1)}, \dots, x_D^{(n)} \in [0, 1]^D$ be the $n$ two-by-two distinct observation points, and $y_1, \dots, y_n \in \mathbb{R}$ be the corresponding observations. Typically, $y_i = f(x_D^{(i)})$ where $f: [0, 1]^D \to \mathbb{R}$ is the function of interest that is modeled by a GP realization. We write $X_D = (x_D^{(1)}, \dots, x_D^{(n)})$ and $y^{(n)} = (y_1, \dots, y_n)$. For $i = 1, \dots, n$ and $\mathcal{J} \subseteq \{1, \dots, D\}$ we write $x_{\mathcal{J}}^{(i)}$ for the vector extracted from $x_D^{(i)}$ by keeping the components with indices in $\mathcal{J}$. We write $X_{\mathcal{J}} = (x_{\mathcal{J}}^{(1)}, \dots, x_{\mathcal{J}}^{(n)})$.

The sequential procedure we suggest is based on updating the set of active input variables and the (multidimensional) subdivision on the space of active input variables. It is initialized by a non-empty set $\mathcal{J}_0 \subseteq \{1, \dots, D\}$ of active variables, with $|\mathcal{J}_0| = d_0$ and by an initial subdivision $S^{(0)} \in \mathcal{S}_{\mathcal{J}_0}$.

We also assume that $x_{\mathcal{J}_0}^{(1)}, \dots, x_{\mathcal{J}_0}^{(n)}$ are two-by-two distinct. This condition is not restrictive and typically holds when selecting a few active variables for initialization (for instance the most correlated with the outputs, empirically). We also assume that the set $\{\alpha \in A_{S^{(0)}}; Y_{S^{(0)},\alpha} \in \mathcal{I}_{X_{\mathcal{J}_0},y^{(n)}} \cap \mathcal{C}_{\mathcal{J}_0}\}$ is non-empty, which means that there exists a finite-dimensional interpolating function satisfying the inequality constraints at the initialization of the sequential procedure. This holds for general choices of a sufficient number of multidimensional knots.

Then, $\hat{\alpha}_{\mathcal{J}_0,S^{(0)},X_{\mathcal{J}_0},y^{(n)}}$ is well-defined when Condition 1 holds, since the minimization set in (8) is non-empty, as discussed above.

### 3.2 The $L^2$ difference between modes

For a non-empty set $\mathcal{J} \subseteq \{1, \dots, D\}$, and for $S \in \mathcal{S}_J$, define

$$ \hat{\alpha}_{\mathcal{J},S} = \hat{\alpha}_{\mathcal{J},S,X_{\mathcal{J}},y^{(n)}} \quad \text{and} \quad \hat{Y}_{\mathcal{J},S} = Y_{S,\hat{\alpha}_{\mathcal{J},S}}. $$
By definition of $\hat{\alpha}_{\mathcal{J},S}$, the function $\hat{Y}_{\mathcal{J},S}$ satisfies both interpolation and inequality constraints, i.e. belongs to $\mathcal{I}_{X_{\mathcal{J},y^{(n)}} \cap \mathcal{C}_{\mathcal{J}}}$. Recall that, by definition of $\hat{\alpha}_{\mathcal{J},S}$, this function is the most likely constrained interpolator of the observations, given the active variables $\mathcal{J}$ and the knots in $S$. Let us now introduce two simple operations on a given subdivision $S \in S_{\mathcal{J}}$.

* Insertion of a new knot: for $i \in \mathcal{J}$ and for $t \in [0, 1] \setminus S_i$, we denote by $S \cup_i t$ the subdivision $S' \in S_{\mathcal{J}}$ defined by $S'_j = S_j$ for $j \in \mathcal{J} \setminus \{i\}$, and by $S'_i = S_i \cup \{t\}$.

* Addition of a new variable: for $i \notin \mathcal{J}$, we denote by $S + i$ the subdivision $S' \in S_{\mathcal{J} \cup \{i\}}$ defined by $S'_j = S_j$ for $j \in \mathcal{J}$, and by $S'_i = S^0$, the minimal subdivision defined in Section 2 that corresponds to functions that are linear with respect to variable $i$.

For a non-empty $\mathcal{J} \subset \{1, \dots, D\}$, for $S \in S_{\mathcal{J}}$ and for $i \in \{1, \dots, D\}$, let

$$N_{S, \mathcal{J}, i} = \begin{cases} \prod_{j \in \mathcal{J} \setminus \{i\}} (|S_j| - 2) & \text{if } i \in \mathcal{J} \\ \prod_{j \in \mathcal{J}} (|S_j| - 2) & \text{if } i \notin \mathcal{J}, \end{cases}$$

with the convention $\prod_{j \neq i} (|S_j| - 2) = 1$ if $\mathcal{J} = \{i\}$.

The quantity $N_{S,\mathcal{J},i}$ is the increase of the number of basis functions if, starting from a subdivision $S$ and a set of active variables $\mathcal{J}$, a knot was inserted for the active variable $i$ or the inactive variable $i$ was made active. Recall that for $j \in \mathcal{J}$, $|S_j| - 2$ is the number of one-dimensional basis functions corresponding to the variable $j$.

Now, we introduce the criterion used by MaxMod. The idea is to measure the $L^2$ difference between the mode of the finite-dimensional process obtained at a potential next step of the algorithm, and the current one, characterized by $\mathcal{J}$ and $S$. At a new iteration, two choices are possible. Either a new knot $t$ is inserted for an active variable $i \in \mathcal{J}$, or a new active variable $i \notin \mathcal{J}$ is added. In this case, the minimal subdivision is added for the dimension $i$ (corresponding to a linear function with respect to the dimension $i$). Formally, for $t \in [0, 1]$ and $i \in \{1, \dots, D\}$, the criterion is written:

$$I_{\mathcal{J},S}(i,t) = \begin{cases} \frac{1}{N_{S,\mathcal{J},i}} \int_{[0,1]^d} (\hat{Y}_{\mathcal{J},S \cup_i t}(x) - \hat{Y}_{\mathcal{J},S}(x))^2 dx & \text{if } i \in \mathcal{J}, \\ \frac{1}{N_{S,\mathcal{J},i}} \int_{[0,1]^{d+1}} (\hat{Y}_{\mathcal{J} \cup \{i\},S+i}(x) - \hat{Y}_{\mathcal{J},S}(x))^2 dx & \text{if } i \notin \mathcal{J}. \end{cases} \quad (10)$$

Note that in (10), in the case $i \notin \mathcal{J}$, we have made the slight abuse of notation of treating $\hat{Y}_{\mathcal{J},S}$ as a function of the $|\mathcal{J}| + 1$ variables $(x_j)_{j \in \mathcal{J} \cup \{i\}}$, that does not use the variable $x_i$. We also remark that, for $i \notin \mathcal{J}$, $I_{\mathcal{J},S}(i,t)$ does not depend on $t$. We explain how to compute (10) efficiently in practice in Appendix B, resulting in a computational complexity that is linear in the number of multi-dimensional knots.

The criterion (10) also penalizes insertions of knots at active variables or additions of new active variables, that increase the number of basis functions (that is the computational complexity) significantly. More precisely, the $L^2$ difference between modes is divided by the number of additional basis functions.

### 3.3 The MaxMod algorithm

For $\Delta, \Delta' > 0$, we introduce the reward

$$R_{\mathcal{J},S}(i,t) = \begin{cases} \Delta d(t, S_i) & \text{if } i \in \mathcal{J} \\ \Delta' & \text{otherwise} \end{cases} \quad (11)$$

that promotes the addition of new variables or the insertion of one-dimensional knots not too close to existing ones. Here, for $q \in \mathbb{N}$, $x \in \mathbb{R}^q$ and $B \subseteq \mathbb{R}^q$, we denote by $d(x, B) = \inf_{u \in B} ||x - u||$, the distance between $x$ and the set $B$. The reward for adding a new active variable is $\Delta'$ and the reward for inserting a knot for an existing variable is $\Delta$ times the distance to the closest existing knot.

For technical reasons, we need to prevent inserting one-dimensional knots that coincide with existing ones. We thus introduce a sequence of separation distances $(b_m)_{m \in \mathbb{N}}$, with $b_m > 0$, such that, at step $m$ of MaxMod, no one-dimensional knot should be inserted at distance less than $b_m$ to an existing knot. In practice, $b_m$ can be taken as small as desired, even equal to the machine precision. We assume that the $(b_m)_{m \in \mathbb{N}}$ are small enough such that

$$\forall m \in \mathbb{N}, \quad 2b_m \left( \max_{i \in \mathcal{J}_0} |S_i^{(0)}| + m \right) < 1. \quad (12)$$
**Algorithm 1** MaxMod (maximum modification of the MAP)

**Input parameters:** $\Delta > 0$, $\Delta' > 0$, two sequence of strictly positive numbers $(a_m)_{m \in \mathbb{N}}$ and $(b_m)_{m \in \mathbb{N}}$, the initial set of active variables $\mathcal{J}_0 \subseteq \{1, \dots, D\}$ and the initial subdivision $S^{(0)} \in S_{\mathcal{J}_0}$.

**Sequential procedure:** For $m \in \mathbb{N}$, $m \ge 0$, do the following.

1: Set $i_{m+1}^* \in \{1, \dots, D\}$, $t_{m+1}^* \in [0, 1]$ such that $d(t_{m+1}^*, S_{i_{m+1}^*}^{(m)}) \ge b_m$ if $i_{m+1}^* \in \mathcal{J}_m$, and such that

$$ 
\begin{aligned}
& I_{\mathcal{J}_m, S^{(m)}}(i_{m+1}^*, t_{m+1}^*) + R_{\mathcal{J}_m, S^{(m)}}(i_{m+1}^*, t_{m+1}^*) + a_m \\
& \ge \sup_{\substack{i \in \{1, \dots, D\}, t \in [0, 1], \\ \text{s.t. } d(t, S_i^{(m)}) \ge b_m \\ \text{if } i \in \mathcal{J}_m}} (I_{\mathcal{J}_m, S^{(m)}}(i, t) + R_{\mathcal{J}_m, S^{(m)}}(i, t)).
\end{aligned} 
\quad (13) $$

2: **if** $i_{m+1}^* \in \mathcal{J}_m$ **then** $\mathcal{J}_{m+1} = \mathcal{J}_m$ and $S^{(m+1)} = S^{(m)} \cup_{i_{m+1}^*} t_{m+1}^*$, that is we add the knot $t_{m+1}^*$ to the current subdivision $S^{(m)}$ for the variable $i_{m+1}^*$.

3: **else** $\mathcal{J}_{m+1} = \mathcal{J}_m \cup \{i_{m+1}^*\}$ and $S^{(m+1)} = S^{(m)} + i_{m+1}^*$, that is we add the variable $i_{m+1}^*$ to the current subdivision $S^{(m)}$.

As will be explained below, (12) simply guarantees that the separation distance does not prevent MaxMod to insert knots.

We also introduce a sequence $(a_m)_{m \in \mathbb{N}}$ of strictly positive numbers that correspond to distances to optimality. The principle is that a knot or a variable added by MaxMod at step $m$ should maximize a quality criterion and that we allow for a distance to the global maximum that is not exactly zero but is simply bounded by $a_m$. The MaxMod sequential procedure can now be written as in Algorithm 1.

MaxMod maximizes a quality criterion, over the added variable or the inserted one-dimensional knot. This quality criterion is the sum of the reward and the $L^2$ distance between the current mode function and the next one (divided by the number of additional knots). Note that if a new variable is added, there is no new one-dimensional knot to select. Indeed, the new one-dimensional subdivision of the new variable is always composed of the knots $\{-1, 0, 1, 2\}$. If a knot is inserted to an existing variable, its location in $[0, 1]$ is optimized continuously. As discussed above, we allow for an approximate maximization of the quality criterion with a gap $a_m$ to the maximum at step $m$. Again as discussed above, no knot should be inserted to an active variable at distance less than $b_m$ from an already existing knot.

Notice that Algorithm 1 is well-defined in the sense that $i_{m+1}^*, t_{m+1}^*$ can be chosen at each step. Indeed, first the supremum in (13) is over a non-empty set. This is because (12) implies that for all $m \in \mathbb{N}$, one can take any $i \in \mathcal{J}_0 \subseteq \mathcal{J}_m$ and find $t \in [0, 1]$ s.t. $d(t, S_i^{(m)}) \ge b_m$. Indeed, consider the intervals of length $2b_m$ centered at $t_j^{(S_i^{(m)})}$, $j = 1, \dots, m$, $S_i^{(m)}$. At step $m$ of Algorithm 1, the number of these intervals is less than $\max_{i \in \mathcal{J}_0} |S_i^{(0)}| + m$, since the initial number of knots is less than $\max_{i \in \mathcal{J}_0} |S_i^{(0)}|$ and the algorithm inserted at most one knot for coordinate $i$ at each previous step. Thus, the union of these intervals does not cover $[0, 1]$. Secondly, the sup is finite in (13) as proved in the following lemma.

**Lemma 1.** Let Condition 1 hold. Then, for each $m \in \mathbb{N}$, the sup in (13) is finite.

# 4 Convergence results

In Section 4.1, we first provide an intermediary result that may be considered of independent interest. We show that for any sequence of multi-dimensional knots, the MAP function converges uniformly to a limit function which we define. Then, in Section 4.2 we use this intermediary result to prove the convergence of MaxMod.

## 4.1 Convergence of the finite-dimensional MAP function with a general sequence of multi-dimensional knots

In this section, we consider a fixed set of active variables $\mathcal{J}$ with $|\mathcal{J}| = d$, and a corresponding sequence of subdivisions, not necessarily obtained from MaxMod. Without loss of generality, we
<table><caption>Table 1: List of symbols for Sections 2 and 3, with the page numbers where the symbols are introduced.</caption><thead><tr><th>Symbol</th><th>Description</th><th>Page</th></tr></thead><tbody><tr><td>φ<sub>u,v,w</sub></td><td>One-dimensional hat basis function with knots u, v, w</td><td>3</td></tr><tr><td>t<sup>(S)</sup><sub>(0)</sub>, ..., t<sup>(S)</sup><sub>(m<sub>S</sub>+1)</sub></td><td>Ordered knots of a one-dimensional subdivision S</td><td>3</td></tr><tr><td>D, d</td><td>Input space dimension and number of active variables</td><td>3</td></tr><tr><td>J = (a<sub>1</sub>, ..., a<sub>d</sub>)</td><td>Set of active variables</td><td>3</td></tr><tr><td>S = (S<sub>a<sub>1</sub></sub>, ..., S<sub>a<sub>d</sub></sub>)</td><td>(Multi-dimensional) subdivision (vector of one-dim. subdivisions)</td><td>3</td></tr><tr><td>l ∈ L<sub>S</sub></td><td>Multi-index for a subdivision S (indexation of the knots)</td><td>3</td></tr><tr><td>α ∈ A<sub>S</sub></td><td>Sequence of coefficients for a subdivision S (values at knots)</td><td>4</td></tr><tr><td>t<sup>(S)</sup><sub>l</sub></td><td>(d-dimensional) knot for a subdivision S</td><td>4</td></tr><tr><td>φ<sup>(S)</sup><sub>l</sub></td><td>Basis function for a subdivision S (tensorized hat basis functions)</td><td>4</td></tr><tr><td>E<sub>S</sub></td><td>Linear space of the componentwise, piecewise linear functions (equivalently multivariate splines of degree 1) with knots in S</td><td>4</td></tr><tr><td>Y<sub>S,α</sub></td><td>Function in E<sub>S</sub> with values at knots α</td><td>4</td></tr><tr><td>π<sub>S</sub>(f)</td><td>Function in E<sub>S</sub> interpolating f at knots in S</td><td>4</td></tr><tr><td>I<sub>U,v<sup>(n)</sup></sub></td><td>Set of d-dimensional functions whose values at inputs U are v<sup>(n)</sup></td><td>5</td></tr><tr><td>C<sub>D</sub></td><td>Set of D-dimensional functions satisfying the inequality constraints</td><td>5</td></tr><tr><td>C<sub>J</sub></td><td>Set of d-dimensional functions satisfying the inequality constraints, when seen as functions of D variables</td><td>5</td></tr><tr><td>k<sub>D</sub>, k<sub>J</sub></td><td>Covariance functions with input dimension D and d</td><td>5</td></tr><tr><td>ξ<sub>J</sub></td><td>d-dimensional GP with covariance function k<sub>J</sub></td><td>5</td></tr><tr><td>k<sub>J</sub>(S, S)</td><td>Covariance matrix of the values of ξ<sub>J</sub> at the knots in S</td><td>6</td></tr><tr><td>M(C<sub>J</sub>), v(C<sub>J</sub>)</td><td>Matrix and vector for the inequality constraints</td><td>6</td></tr><tr><td>α̂<sub>J,S,U,v<sup>(n)</sup></sub></td><td>Values at knots of ξ<sub>J</sub> with highest density conditionally on interpolation and inequality constraints</td><td>7</td></tr><tr><td>Y<sub>S,α̂</sub><sub>J,S,U,v<sup>(n)</sup></sub></td><td>d-dim. function in E<sub>S</sub> built from α̂<sub>J,S,U,v<sup>(n)</sup></sub> (MAP or mode)</td><td>7</td></tr><tr><td>X<sub>D</sub>, X<sub>J</sub>, y<sup>(n)</sup></td><td>Observation points and associated observed values for MaxMod</td><td>7</td></tr><tr><td>α̂<sub>J,S</sub>, Ŷ<sub>J,S</sub></td><td>α̂<sub>J,S,U,v<sup>(n)</sup></sub> and Y<sub>S,α̂</sub><sub>J,S,U,v<sup>(n)</sup></sub> for MaxMod</td><td>7</td></tr><tr><td>S ∪<sub>i</sub> t</td><td>Insertion of the knot t to the subdivision S for the variable i</td><td>8</td></tr><tr><td>S + i</td><td>Addition of the variable i to the subdivision S</td><td>8</td></tr><tr><td>N<sub>S,J,i</sub></td><td>Number of additional basis functions for a new knot/variable</td><td>8</td></tr><tr><td>I<sub>J,S</sub>(i, t)</td><td>Normalized L<sup>2</sup> difference between modes for a new knot/variable</td><td>8</td></tr><tr><td>R<sub>J,S</sub>(i, t)</td><td>Reward promoting new variables and distance between knots</td><td>8</td></tr><tr><td>d(x, B)</td><td>Distance between a point x and a set B</td><td>8</td></tr></tbody></table>

set $J = \{1, \dots, d\} \subseteq \{1, \dots, D\}$ for $d \le D$. Thus we shall consider functions from $[0, 1]^d$ to $\mathbb{R}$, with the fixed dimension $d$. This allows us to remove the dependence on $J$ in the notations, and to write more simply: $X = (x^{(1)}, \dots, x^{(n)}) = (x_J^{(1)}, \dots, x_J^{(n)})$, $S = S_J$, $\hat{\alpha}_{S,U,v^{(n)}} = \hat{\alpha}_{J,S,U,v^{(n)}}$, $\hat{\alpha}_S = \hat{\alpha}_{J,S}$, $\hat{Y}_S = \hat{Y}_{J,S}$, $C = C_J$ and $k = k_J$. In particular, the sets $C$ corresponding to (3), (4) and (5) are written:

$$ C = \{f \in C([0, 1]^d, \mathbb{R}); a \le f(x) \le b \text{ for all } x \in [0, 1]^d\}, \quad (14) $$

$$ C = \{f \in C([0, 1]^d, \mathbb{R}); f(u) \le f(v) \text{ for all } u, v \in [0, 1]^d, u \le v\} \quad (15) $$

$$ C = \{f \in C([0, 1]^d, \mathbb{R}); \text{ for all } i \in \{1, \dots, d\}, \text{ for all } x_{\sim i} \in [0, 1]^{d-1}, \quad (16) $$

the function $u_i \mapsto f(u_i, x_{\sim i})$ is convex\}.

### 4.1.1 The multiaffine extension of a multivariate function

For a univariate continuous function $g$ defined on a closed subset $B$ of $[0, 1]$ containing 0 and 1, let us denote by $L_B(g)$, or simply $L_B g$, the *affine extension* of $g$ on $[0, 1]$:

$$ L_B g(t) = \begin{cases} g(t) & \text{if } t \in B \\ \omega_-(t)g(t^-) + \omega_+(t)g(t^+) & \text{if } t \notin B \end{cases} \quad (17) $$

where $t^- = \max\{u \in B; u \le t\}$ is the closest left neighbor of $t$ in $B$, $t^+ = \min\{u \in B; u \ge t\}$ is the closest right neighbor of $t$ in $B$, $\omega_+(t) = \frac{t-t^-}{t^+-t^-}$ and $\omega_-(t) = 1 - \omega_+(t)$. Notice that, thanks to the
<figure><img src="image_3.png" /><figcaption>Figure 2: Sequential construction of the multiaffine extension. Illustration on the 2-dimensional function $f(x) = (x_1 - 0.5)^2(x_2 - 0.5)^3$, defined on $F_1 \times F_2$, with $F_1 = F_2 = [0, 0.1] \cup [0.4, 1]$. The three perspective plots represent, from left to right: $f$, $L_2f$ and $P_{F \to [0,1]^d}f = L_1L_2f$.</figcaption></figure>

assumptions on $B$, if $t \in [0, 1] \setminus B$ then $t^-$ and $t^+$ are well-defined and distinct and belong to $B$. Hence $\omega_-(t), \omega_+(t) \in (0, 1)$ are well-defined, and $L_B$ is well-defined by (17).

One can show that $L_B g$ is the unique continuous function $h$ equal to $g$ on $B$ such that $h$ is affine on all intervals of $[0, 1] \setminus B$ (see the proof of Proposition 1). This is clear when $B$ is a finite union of intervals, but it is still valid for more complex closed sets, such as an infinite sequence with an accumulation point. Starting from this property, we now introduce the multiaffine extension of a continuous multivariate function, as defined below.

**Definition 1 (multiaffine function).** Let $f$ be a function defined on a product set $G = G_1 \times \dots \times G_d$. We say that $f$ is $d$-affine or simply multiaffine if it is componentwise affine: for all $i = 1, \dots, d$, for all $u = (u_1, \dots, u_d) \in G$, the function $u_i \mapsto f(u)$ is affine on $G_i$.

**Proposition 1 (definition of the multiaffine extension).** Consider a continuous multivariate function $f$ on $F = F_1 \times \dots \times F_d$, where each $F_i$ is a closed subset of $[0, 1]$ containing 0 and 1. Then, there exists a unique continuous function $g$ defined on $[0, 1]^d$, equal to $f$ on $F$, such that, for all $i = 1, \dots, d$, for all $t \in [0, 1]^d$, the univariate cut function $g(\cdot, t_{\sim i}) : u_i \mapsto g(u_i, t_{\sim i})$ is affine on each interval of $[0, 1] \setminus F_i$. Furthermore, $g$ is obtained sequentially from $f$ by linearly extending univariate cuts. More precisely, for $i \in \{1, \dots, d\}$, denote $F^{1:i} = F_1 \times \dots \times F_i$, $F^{1:0} = \emptyset$, and let $L_i$ be defined from $\mathcal{F}(F^{1:i} \times [0, 1]^{d-i}, \mathbb{R})$ to $\mathcal{F}(F^{1:(i-1)} \times [0, 1]^{d-i+1}, \mathbb{R})$ by

$$L_i(h)(t) = L_i h(t) = (L_{F_i}[h(\cdot, t_{\sim i})])(t_i)$$

for $h \in \mathcal{F}(F^{1:i} \times [0, 1]^{d-i}, \mathbb{R})$ and $t \in F^{1:(i-1)} \times [0, 1]^{d-i+1}$. Then

$$g = L_1 \dots L_d f.$$

We call $g$ the multiaffine extension of $f$, and denote it by $P_{F \to [0,1]^d}(f)$ (or simply $P_{F \to [0,1]^d}f$). Furthermore, $P_{F \to [0,1]^d}(f)$ is given explicitly by, for $t \in [0, 1]^d$,

$$P_{F \to [0,1]^d}(f)(t) = \sum_{\epsilon_1, \dots, \epsilon_d \in \{-, +\}} \left( \prod_{j=1}^{d} \omega_{\epsilon_j}(t_j) \right) f(t_1^{\epsilon_1}, \dots, t_d^{\epsilon_d}), \quad (18)$$

where the $2^d$ products are non-negative and sum to one.

A two-dimensional illustration of the multiaffine extension is provided in Figure 2. As observed on this figure, a feature of the multiaffine extension is the extension of the definition of a continuous function $f$ on hyper-rectangles where $f$ is defined only at the $2^d$ vertices (for instance the square $[0.1, 0.4]^2$ in the figure). The connection with the hat basis functions of Section 2 is that the extension at one of these hyper-rectangles coincides with the expression of $Y_{S,\alpha}$ from (1), when the closest knots in $S$ to the hyper-rectangle coincide with its vertices. We refer to the last item of Remark 1 for a more formal statement.
**Remark 1.**

* As a direct consequence of Proposition 1, all $q$-dimensional cuts of the multiaffine extension, i.e. $(t_{i_1}, \dots, t_{i_q}) \mapsto P_{F \to [0,1]^d} f(t)$, are $q$-affine on any hypercube in the product of complementary sets $\prod_{l=1}^q ([0, 1] \setminus F_{i_l})$.

* In Proposition 1, the multiaffine extension is obtained by composing 1-dimensional cuts in a specific order. By uniqueness, permuting the order will give the same result.

* In one dimension, extending linearly a function by using the closest neighbors as in (17) can be viewed as computing the posterior mean of the Brownian motion [10, preface, page xiv]. In general, one could define the multiaffine extension of Proposition 1 with the Brownian sheet. However, this would require to condition the Brownian sheet on a continuous set and would add technicalities (see, e.g., [4] for more details). This is why we have chosen to introduce the multiaffine extension with basic tools.

* Consider $S \in \mathcal{S}$ and the set $F_S = \prod_{j=1}^d (S_j \cap [0, 1])$, where each $S_j \cap [0, 1]$ is a closed subset of $[0, 1]$ containing 0 and 1. Then, by Proposition 2, for all $f \in \mathcal{C}([0, 1]^d, \mathbb{R})$,

$$ \pi_S(f) = P_{F_S \to [0,1]^d}(f|_{F_S}). \qquad (19) $$

Hence, we can also interpret the multiaffine extension as a generalization of the projection $\pi_S$ to an infinite number of knots.

Now, we gather below some properties of the multiaffine extension.

**Proposition 2.**

1. The map $f \mapsto P_{F \to [0,1]^d}(f)$ is affine and 1-Lipschitz from $\mathcal{C}(F, \mathbb{R})$ to $\mathcal{C}([0, 1]^d, \mathbb{R})$, equipped with the $L^\infty$ norm. In particular, it preserves uniform convergence.

2. If $S$ is a $d$-dimensional subdivision such that $S \cap [0, 1]^d \subseteq F$, then any piecewise multilinear function constructed from $S$ coincides with the multiaffine extension of its restriction to $F$:

$$ \forall f \in E_S: \quad P_{F \to [0,1]^d}(f|_F) = f. $$

**Corollary 1.** For $x \in [0, 1]^d$, let $f, g$ be $d$-affine on the hypercube $\Delta = \prod_{j=1}^d [x_j^-, x_j^+]$, and coincide on its $2^d$ vertices. Then they are equal on $\Delta$.

### 4.1.2 Technical conditions

We let $m_0 \in \mathbb{N}$ and consider a sequence of subdivisions $(S^{(m)})_{m \ge m_0}$ with $S^{(m)} \in \mathcal{S}$ for $m \ge m_0$ (not necessarily the sequence obtained from MaxMod). Let the subdivisions be nested (i.e. $S_i^{(m)} \subseteq S_i^{(m+1)}$ for $i \in \{1, \dots, d\}$). We will show the convergence of the sequence of mode functions obtained from the sequence of subdivisions $(S^{(m)})_{m \ge m_0}$. In that view, we now introduce a list of technical conditions.

To the set $\mathcal{C}$, that is interpreted as the set of functions satisfying inequality constraints, we associate a set $\dot{\mathcal{C}} \subseteq \mathcal{C}$, that we can choose and that is interpreted as a set of functions satisfying corresponding strict inequality constraints. Let us explicitly show how $\dot{\mathcal{C}}$ is chosen when $\mathcal{C}$ is given by one of (14), (15) or (16).

First, consider that $\mathcal{C}$ is given by (14), with $-\infty < a < b < +\infty$ (the cases where either $a = -\infty$ or $b = +\infty$ are similar). Then we choose $\dot{\mathcal{C}}$ as the set of continuous functions that are strictly between $a$ and $b$. Second, consider that $\mathcal{C}$ is given by (15). Let us write $u < v$ when $u \le v$ and $u \ne v$. We say that a function $f$ from $[0, 1]^d$ to $\mathbb{R}$ is strictly increasing if $f(u) < f(v)$ for all $u, v \in [0, 1]^d$, $u < v$. Then we chose $\dot{\mathcal{C}}$ as the set of continuous strictly increasing functions. Finally, consider that $\mathcal{C}$ is given by (16). Then we choose $\dot{\mathcal{C}}$ as

$$ \{f \in \mathcal{C}([0, 1]^d, \mathbb{R}), \quad \begin{array}{l} \text{for all } i \in \{1, \dots, d\}, \text{ for all } x_{\sim i} \in [0, 1]^{d-1}, \\ \text{the function } u_i \mapsto f(u_i, x_{\sim i}) \text{ is strictly convex.} \end{array} \} $$

For other inequality sets $\mathcal{C}$, we emphasize that we are free to choose the set $\dot{\mathcal{C}} \subseteq \mathcal{C}$ for which the technical conditions given below (Conditions 2 and 4) hold.

**Condition 2** (initial knots flexibility, 1). The set $\{\alpha \in A_{S^{(m_0)}}, Y_{S^{(m_0)}, \alpha} \in \mathcal{I}_{X, y^{(n)} \cap \dot{\mathcal{C}}}\}$ is non-empty.
The previous condition means that one can construct, from the initial subdivision, a finite-
dimensional function that satisfies the interpolation and strict inequality constraints. With the
choice of $\mathcal{C}$ for boundedness, Condition 2 implies that $a < y_1 < b, \dots, a < y_n < b$ and that there
are enough knots in the initial subdivision to generate an interpolating function that is strictly
between $a$ and $b$.

With the choice of $\mathcal{C}$ for monotonicity, Condition 2 implies $y_i < y_j$ for $x_i < x_j$ ($i, j = 1, \dots, n$) and that there are sufficiently many knots in the initial subdivision to generate a strictly increasing interpolating function.

Consider finally componentwise convexity. For $U = (u_1, \dots, u_r) \in ([0, 1]^d)^r$ and $v^{(r)} = (v_1, \dots, v_r) \in \mathbb{R}^r$, let us say that $U$ and $v^{(r)}$ are compatible with strict convexity if there exists a function $g$ in $\mathcal{C}$ such that $g(u_i) = v_i, i = 1, \dots, r$. Then, Condition 2 implies that $X$ and $y^{(n)}$ are compatible with strict convexity and that there are sufficiently many knots in the initial subdivision to generate a componentwise strictly convex interpolating function.

**Condition 3** (initial knots flexibility, 2). $\{ (Y_{S^{(m_0)}, \alpha}(x^{(i)}))_{i=1,\dots,n} ; \alpha \in S^{(m_0)} \} = \mathbb{R}^n$.

This condition means that for any possible $n$-dimensional observation vector (not only $y^{(n)}$), there exists a combination of basis functions based only on the initial knots that interpolates this observation vector. This condition requires that there are enough initial knots (at least $n$) and that these knots are located adequately. This condition is mild, since the number of knots increases to infinity. For instance, in dimension $d = 1$, it is sufficient that there is one initial knot in between each pair of observation points.

For a subset $B$ in a metric space with distance dist, its interior is written int$_{dist}$(B) and its closure is written $\bar{B}$. In particular, for a finite dimensional subspace $A$, and for $B \subseteq C(A, \mathbb{R})$, we shall consider int$_{||\cdot||_\infty}$(B) relatively to the $L^\infty$ norm. We let $\mathcal{H}$ be the RKHS (see for instance [10]) of $k$ and we write $||\cdot||_\mathcal{H}$ for the RKHS norm on $\mathcal{H}$. We shall then consider int$_{||\cdot||_\mathcal{H}}(\mathcal{H} \cap C)$, where the interior is defined w.r.t. the RKHS norm of $k$ on $\mathcal{H}$.

**Condition 4** (RKHS interior). For all $h \in \mathcal{C}$, for all $r \in \mathbb{N}$, $U = (u_1, \dots, u_r) \in ([0, 1]^d)^r$, with $u_1, \dots, u_r$ two-by-two distinct, letting $v^{(r)} = (h(u_1), \dots, h(u_r)) \in \mathbb{R}^r$, the set int$_{||\cdot||_\mathcal{H}}(\mathcal{H} \cap C) \cap \mathcal{I}_{U,v^{(r)}}$ is non-empty.

With the choice of $\mathcal{C}$ for boundedness, Condition 4 holds when int$_{||\cdot||_\mathcal{H}}(\mathcal{H} \cap C) \cap \mathcal{I}_{U,v^{(r)}}$ is non-empty for all $U = (u_1, \dots, u_r) \in ([0, 1]^d)^r$ and $v^{(r)} = (v_1, \dots, v_r) \in \mathbb{R}^r$, with $u_1, \dots, u_r$ two-by-two distinct and $a < v_1 < b, \dots, a < v_n < b$. We also have the following lemma.

**Lemma 2.** With $C$ given by (14) with $-\infty < a < b < +\infty$ and when $\mathcal{C}$ is the set of continuous functions that are strictly between $a$ and $b$, we have $\mathcal{H} \cap \mathcal{C} \subseteq \text{int}_{||\cdot||_\mathcal{H}}(\mathcal{H} \cap C)$.

From Lemma 2, the set int$_{||\cdot||_\mathcal{H}}(\mathcal{H} \cap C) \cap \mathcal{I}_{U,v^{(r)}}$ is non-empty if there exists $h \in \mathcal{H}$ that is strictly between $a$ and $b$ and interpolates $(u_1, v_1), \dots, (u_v, v_n)$. This can be interpreted as requiring $\mathcal{H}$ to be rich enough, and is not restrictive (it is also required in [5, 6]).

With the choice of $\mathcal{C}$ for monotonicity, Condition 4 holds when the set int$_{||\cdot||_\mathcal{H}}(\mathcal{H} \cap C) \cap \mathcal{I}_{U,v^{(r)}}$ is non-empty for all $U = (u_1, \dots, u_r) \in ([0, 1]^d)^r$ and $v^{(r)} = (v_1, \dots, v_r) \in \mathbb{R}^r$, with $u_1, \dots, u_n$ two-by-two distinct, with $v_i < v_j$ for $u_i < u_j$ ($i, j = 1, \dots, n$). Again, this is not restrictive, it is also required in [5, 6] and can be interpreted as $\mathcal{H}$ being rich enough.

With the choice of $\mathcal{C}$ for componentwise convexity, Condition 4 holds when int$_{||\cdot||_\mathcal{H}}(\mathcal{H} \cap C) \cap \mathcal{I}_{U,v^{(r)}}$ is non-empty for all $U = (u_1, \dots, u_r) \in ([0, 1]^d)^r$ and $v^{(r)} = (v_1, \dots, v_r) \in \mathbb{R}^r$ that are compatible with strict convexity. Again this is not restrictive and is also required in [5, 6].

Next, we introduce two conditions related to the stability of the constraint set.

**Condition 5** (constraint set stability by projection). For all subdivision $S \in \mathcal{S}$, $\pi_S(C) \subseteq C$.

**Condition 6** (constraint set stability by multiaffine extension). For $f \in C$, with $f|_F$ the restriction of $f$ to $F$, we have $P_F \to [0,1]^d$ ($f|_F \in C$).

By Lemma 3 below, these two conditions hold in the cases of boundedness, monotonicity and componentwise convexity. About Condition 5, this is a consequence of Remark 1 and Lemma 3, applied to the set $F_S = S \cap [0,1]^d$, that is closed, of product form, and contains $\{0,1\}^d$. Notice that the fact that Condition 5 holds for boundedness and monotonicity constraints has been already proved in [42].

**Lemma 3.** When $C$ is given by one of (14), (15) or (16), Condition 6 holds.
### 4.1.3 Convergence of the finite-dimensional MAP function

For a sequence of sets $(B_m)_{m \ge m_0}$, with $B_m \subseteq [0, 1]$ for $m \ge m_0$, we say that $B_m$ is dense in $[0, 1]$ as $m \to \infty$ if for every $x \in [0, 1]$, we have $d(x, B_m) \to 0$ as $m \to \infty$. We first recall the convergence theorem (Theorem 3) in [6] that shows that, when for $i = 1, \dots, d$ the knots in $S_i^{(m)}$ are dense in $[0, 1]$, then the function $\hat{Y}_{S_i^{(m)}}$ converges as $m \to \infty$ to a limit function $Y_{opt}$.

**Theorem 1** (Kimeldorf-Wahba correspondence under constraints, [6]). Consider a sequence of nested subdivisions $(S^{(m)})_{m \ge m_0}$ with $S^{(m)} \in S$ for $m \ge m_0$ (i.e. $S_i^{(m)} \subseteq S_i^{(m+1)}$ for $i \in \{1, \dots, d\}$), such that for $i \in \{1, \dots, d\}$, $S_i^{(m)}$ is dense in $[0, 1]$. Assume that Condition 1 holds. Assume that the two following conditions hold.

$(H_1)$ $\text{int}_{||\cdot||_H}(\mathcal{H} \cap C) \cap \mathcal{I}_{X,y^{(n)}}$ is non-empty,

$(H_2)$ $\forall m \ge m_0, \pi_{S^{(m)}}(C) \subseteq C$.

Then, as $m \to \infty$, the function $\hat{Y}_{S^{(m)}}$ converges uniformly on $[0, 1]^d$ to $Y_{opt}$, with:

$$Y_{opt} = \underset{f \in \mathcal{H} \cap \mathcal{C} \cap \mathcal{I}_{X,y^{(n)}}}{\operatorname{argmin}} \|f\|_{\mathcal{H}}.$$

In fact, in the above theorem, $(H_1)$ and $(H_2)$ hold in our framework. Indeed, $(H_1)$ is a consequence of Conditions 2 and 4 and $(H_2)$ is a consequence of Condition 5.

The interpretation is that $\hat{Y}_{S^{(m)}}$ is a finite-dimensional GP mode, obtained by solving an optimization problem in finite dimension (see Section 2), while $Y_{opt}$ is the optimal constrained interpolator in $\mathcal{H}$, obtained by solving an optimization problem in infinite dimension, see [6].

We now give an extension to Theorem 1 to the case where, for $i \in \{1, \dots, d\}$, the knots in $S_i^{(m)}$ are not necessarily dense in $[0, 1]$. This extension is the main result of this section. We define $F = F_1 \times \dots \times F_d$, where $F_i$ is the closure in $[0, 1]$ of all the knots at coordinate $i$:

$$F_i = [0, 1] \cap \overline{\bigcup_{m \ge m_0} S_i^{(m)}}.$$

We further denote by $\mathcal{C}_F$ and $\mathcal{I}_{F,X,y^{(n)}}$ the pre-image sets of $\mathcal{C}$ and $\mathcal{I}_{X,y^{(n)}}$ by the multiaffine extension $P_{F \to [0,1]^d}$ introduced in Section 4.1.1:

$$\begin{aligned} \mathcal{C}_F &= \{f : F \to \mathbb{R}, \text{ s.t. } P_{F \to [0,1]^d} f \in \mathcal{C}\}, \\ \mathcal{I}_{F,X,y^{(n)}} &= \{f : F \to \mathbb{R}, \text{ s.t. } P_{F \to [0,1]^d} f \in \mathcal{I}_{X,y^{(n)}}\}. \end{aligned}$$

We also let $k_F$ be the restriction of the kernel $k$ on $F \times F$ and $\mathcal{H}_F$ be the corresponding RKHS of functions from $F \to \mathbb{R}$ (see e.g. [10]). We have $\mathcal{H}_F = \{f : F \to \mathbb{R}, \exists h \in \mathcal{H} \text{ s.t. } h|_F = f\}$, and the RKHS norm in $\mathcal{H}_F$ is $\|f\|_{\mathcal{H}_F} = \inf_{h|_F=f} \|h\|_{\mathcal{H}}$.

Then we state the extension of Theorem 1 to non-dense sequences. The result is intuitive: on the closure set $F$, the finite-dimensional GP mode converges uniformly to the constrained interpolator in $\mathcal{H}_F$, with the equality and inequality constraints given by $\mathcal{I}_{F,X,y^{(n)}}$ and $\mathcal{C}_F$. On the complement of $F$, these two functions are piecewise multilinear (the interpolator in $\mathcal{H}_F$ being extended with the multiaffine extension), see Proposition 1. The multiaffine extension enables to express the convergence, both on $F$ and its complement, simply.

**Theorem 2.** Consider a sequence of nested subdivisions $S^{(m)}$, as well as $F$, $\mathcal{H}_F$, $\mathcal{C}_F$, $\mathcal{I}_{F,X,y^{(n)}}$, as defined in this section. Assume that Conditions 1 to 6 hold. Then, as $m \to \infty$, the function $\hat{Y}_{S^{(m)}}$ converges uniformly on $[0, 1]^d$ to $P_{F \to [0,1]^d}(Y_{F,\text{opt}})$, with:

$$Y_{F,\text{opt}} = \underset{f \in \mathcal{H}_F \cap \mathcal{C}_F \cap \mathcal{I}_{F,X,y^{(n)}}}{\operatorname{argmin}} \|f\|_{\mathcal{H}_F}.$$

Note that the two functions $Y_{opt}$ in Theorem 1 and $Y_{F,\text{opt}}$ in Theorem 2 are equal when $F = [0, 1]^d$. When $F \neq [0, 1]^d$, these two functions need not coincide, even on $F$.

## 4.2 Convergence of the MaxMod algorithm

We can now apply Theorem 2 to prove the convergence of MaxMod in Theorem 3 below. The technical conditions for Theorem 3 are adaptations of those of Theorem 2 to the setting of MaxMod.
They are stated and discussed in Appendix C. There, we also show that these conditions hold in the cases of boundedness, monotonicity and componentwise convexity.

Theorem 3 shows the consistency of MaxMod, which will asymptotically select all the variables and allocate a dense sequence of knots to each variable. As a consequence, the mode function obtained from MaxMod converges to the infinite-dimensional optimal function $Y_{opt}$ defined in Theorem 1. In Theorem 3, we let $\mathcal{H}_D$ be the RKHS of $k_D$.

**Theorem 3.** Let $(\mathcal{J}_m)$ and $(S^{(m)})$ be the sequence of sets of active variables and of subdivisions obtained from Algorithm 1. Assume that $a_m \to 0$ and that Conditions 1 and 7 to 11 (see Appendix C) hold. Then, for $m$ large enough, $\mathcal{J}_m = \{1, \dots, D\}$. Furthermore, for $j = 1, \dots, D$, the set $S_j^{(m)}$ (which becomes well-defined for $m$ large enough) is dense in $[0, 1]$ as $m \to \infty$. Consequently, as $m \to \infty$, the mode $\hat{Y}_{S^{(m)}}$ converges uniformly on $[0, 1]^D$ to the function $Y_{opt}$ defined by

$$Y_{opt} = \underset{f \in \mathcal{H}_D \cap \mathcal{C}_D \cap \mathcal{I}_{X_D, y^{(n)}}}{\operatorname{argmin}} \|f\|_{\mathcal{H}_D}.$$

Theorem 3 can be interpreted as stating that, for a fixed dataset, as the computational budget (quantified here by the number of multi-dimensional knots) goes to infinity, the finite-dimensional mode converges to $Y_{opt}$, which is optimal for the dataset, but requires, so to speak, an infinite computational budget. Remark that MaxMod is sequential and of the greedy type. Hence, as discussed in Section 1, it is important to guarantee that its one-step-ahead allocation of the knots does not prevent it to yield mode functions that are converging to the global optimum $Y_{opt}$.

**Remark 2.** From the proof of Theorem 3, one can see that the convergence of MaxMod still holds if the $L^2$ distance in (10) is replaced by any discrepancy criterion $\Delta(f_1, f_2)$, that goes to zero when $\|f_1 - f_2\|_{\infty} \to 0$.

# 5 Numerical experiments

In this section, we aim at testing the performance of the constrained GP when the knots and active dimensions are sequentially added using MaxMod. In practice, as shown in [38], incorporating noise in the constrained GP model leads to significant computational improvements due to the “relaxation” of the interpolation conditions. The noise, parametrized by a variance $\tau^2 \ge 0$, leads to a new definition of the mode in (9), that can be found in [38, 41], and that we call the noisy mode. This noisy mode can be computed even when the number of multi-dimensional knots $m = m_{S_{a_1}} \times \dots \times m_{S_{a_d}}$ is smaller than the number of observations $n$. MaxMod can be carried out exactly as in Section 3 when working with the noisy mode. We will denote as $\hat{Y}_{MaxMod}$ the noisy mode obtained from MaxMod.

Here, we shall work with the noisy mode, which allows us to always initialize MaxMod with only one active dimension, i.e. $|\mathcal{J}| = 1$, and to add new active ones according to Algorithm 1. We take the first active dimension as the one resulting in an initial mode $\hat{Y}_{MaxMod,0}$ that differs the most from zero. Then, for each addition of a new dimension, an initial set of two knots is allocated at the boundaries of $[0, 1]$.

Algorithm 2 summarizes the practical implementation of MaxMod that is used within this section. We fix $\Delta = \Delta' = 1 \times 10^{-9}$, see (11), which means that there is a negligible reward for either inserting a new knot in an already active dimension or adding a new dimension. For other applications, one may be particularly interested in adding knots while preserving a tractable input dimension, in which case $\Delta$ must be larger than $\Delta'$. Furthermore, we consider a squared exponential kernel for the covariance function $k_{\mathcal{J}}$ [23, 45]. For each step of MaxMod, the covariance parameters of $k_{\mathcal{J}}$ and the noise variance parameter $\tau^2$ are estimated via maximum likelihood [45]. The open source code of MaxMod is available in the R package lineqGPR (0.2.0) [37].<fnref n="1" />

## 5.1 2D illustration under monotonicity constraints

For illustration, we consider the 2D monotonic function $f(x) = \frac{1}{2}x_1 + \arctan(10x_2)$ on $[0, 1]^2$. We evaluate $f$ at a maximin Latin hypercube design (LHD, [19]) with $n = 40$ points $x_D^{(1)}, \dots, x_D^{(n)}$.

<footnote n="1">¹The source code can be found at: https://github.com/anfelopera/lineqGPR. The analytic examples proposed here (Sections 5.1 and 5.2) can be reproduced using the Jupyter notebooks available in the GitHub repository (https://github.com/anfelopera/lineqGPR/tree/master/notebooks). The numerical experiments were executed on an 11th Gen Intel(R) Core(TM) i5-1145G7 CPU@2.60GHz, 16 Gb RAM. Computations from step 2 to 4 of Algorithm 2 were parallelized using two independent cores.</footnote>
**Algorithm 2** Practical implementation of MaxMod using the noisy mode

**Input parameters:** $\Delta > 0$, $\Delta' > 0$, $D$.

1: Initialize the algorithm by selecting as initial input dimension the variable that gives the largest norm of the initial noisy mode $\hat{Y}_{MaxMod,0}$.

**Sequential procedure:** For $m \in \mathbb{N}$, $m \ge 0$, do the following.

2: **for** $k = 1, \dots, D$ **do**

3: **if** the variable $k$ is already active **then** compute the optimal position of the new knot $t_k \in [0, 1]$ according to the MaxMod criterion described in (13) (in the supremum), and denote the resulting noisy mode as $\hat{Y}_{MaxMod,m+1}^{(k)}$.

4: **else** add two knots at the boundaries of the selected new active dimension, i.e. $(t_{k,1}, t_{k,2}) = (0, 1)$, and denote the resulting noisy mode as $\hat{Y}_{MaxMod,m+1}^{(k)}$.

5: Choose the optimal decision $k^* \in \{1, \dots, D\}$ that maximizes the MaxMod criterion:

$$k^* \in \underset{k \in \{1, \dots, D\}}{\operatorname{argmax}} \operatorname{MaxMod}_{\Delta, \Delta'}(\hat{Y}_{\operatorname{MaxMod},m}, \hat{Y}_{\operatorname{MaxMod},m+1}^{(k)}).$$

6: Update knots and active variables, set new noisy mode to $\hat{Y}_{MaxMod,m+1} = \hat{Y}_{MaxMod,m+1}^{(k^*)}$.

<figure><img src="image_4.png" /><figcaption>Figure 3: Evolution of MaxMod using $f(x) = \frac{1}{2}x_1 + \arctan(10x_2)$ as target function. The mode $\hat{Y}_{MaxMod}$ accounts for monotonicity constraints everywhere. The panels shows: the observations (black dots), the mode (1D: green solid line, 2D: solid surface) and the knots (red crosses). The set of added knots are highlighted by a vertical red line in 1D and a vertical red plane in 2D.</figcaption></figure>

The observations, with the notation of Section 3, are thus $y_1 = f(x_D^{(1)}), \dots, y_n = f(x_D^{(n)})$. For $\hat{Y}_{MaxMod}$, we account for monotonicity constraints everywhere. As a stopping rule, we check that the criterion in (10), plus the reward in (11), is smaller than a tolerance fixed to $1 \times 10^{-5}$, for all possible new knot or variable.

### 5.1.1 Assessment of MaxMod in terms of bending energy

From Figure 3, we can observe that MaxMod starts by adding the second dimension rather than the first one since $f$ is more variable across $x_2$. Note also that, before activating the first dimension, the algorithm refines the second one by placing a third knot around $x_2 = 0.29$. Then, after the second iteration, although the first dimension has been activated, the algorithm prefers enriching the quality of $\hat{Y}_{MaxMod}$ across $x_2$ while staying linear across $x_1$. The convergence of the algorithm is obtained after four iterations, resulting in a total of $m = 10$ knots: 2 and 5 one-dimensional knots allocated across the first and the second dimension, respectively. We note that the final estimated noise variance parameter is negligible, $\hat{\tau}_{MaxMod}^2 = 8.38 \times 10^{-5}$ (equivalent to $\approx 0.07\%$ of the variance of the observations), resulting in a GP model that almost interpolates all the observations (see Figure 3).

We now compare the mode $\hat{Y}_{MaxMod}$ to modes resulting from equispaced designs of the knots (see [38, 39, 42]). We consider either square or rectangular designs of the knots. This leads to three modes denoted as:
<figure><img src="image_5.png" /><figcaption>Figure 4: Evolution of the (normalized) bending energy $E_n$ for the example in Figure 3. Results are shown for $\hat{Y}_{\text{square}}$ (green dotted line), $\hat{Y}_{\text{MaxMod,rect}}$ (red dashed line) and $\hat{Y}_{\text{MaxMod}}$ (blue solid line). For $\hat{Y}_{\text{MaxMod,rect}}$ and $\hat{Y}_{\text{MaxMod}}$, the labels on top are given for each iteration. They denote which dimension has been refined by MaxMod and vertical dashed lines indicate when a new active dimension is added.</figcaption></figure>

* $\hat{Y}_{\text{MaxMod}}$: the mode resulting from MaxMod.

* $\hat{Y}_{\text{MaxMod, rect}}$: the mode resulting from equispaced one-dimensional knots but where the number of knots per dimension is the same as for $\hat{Y}_{\text{MaxMod}}$.

* $\hat{Y}_{\text{square}}$: the mode resulting from equispaced one-dimensional knots with the same number of knots per dimension.

We assess the quality of $\hat{Y}_{\text{square}}$, $\hat{Y}_{\text{MaxMod,rect}}$ and $\hat{Y}_{\text{MaxMod}}$ in terms of the (normalized) bending energy [17]: $E_n(f, \hat{Y}) = \int_{[0,1]^D} (f(x) - \hat{Y}(x))^2 dx / \int_{[0,1]^D} f^2(x) dx$. By comparing the $E_n$ values using $\hat{Y}_{\text{MaxMod}}$ to the ones provided by $\hat{Y}_{\text{square}}$ or $\hat{Y}_{\text{MaxMod,rect}}$, we aim at showing that MaxMod not only adds active dimensions strategically but also places knots in regions leading to smaller errors.

Figure 4 shows the performance of $\hat{Y}_{\text{square}}$, $\hat{Y}_{\text{MaxMod,rect}}$ and $\hat{Y}_{\text{MaxMod}}$ in terms of $E_n$. Observe that the mode $\hat{Y}_{\text{MaxMod}}$ minimizes faster the $E_n$ criterion, leading to negligible values after adding $m = 8$ knots (iteration 4, see also Figure 3). For equispaced designs of knots, the $E_n$ values from $\hat{Y}_{\text{MaxMod,rect}}$ outperformed those from $\hat{Y}_{\text{square}}$. However, both $\hat{Y}_{\text{square}}$ and $\hat{Y}_{\text{MaxMod,rect}}$ led to suboptimal results due to the equispaced restriction.

### 5.1.2 Computation time

We report the computation time taken by MaxMod and compare it with equispaced knots evenly allocated to each variable (mode $\hat{Y}_{\text{square}}$). We consider the execution of MaxMod up to $m = 6$ and the equispaced approach with $m = 16$. Indeed, from Figure 4, both models lead to similar $E_n$ values. MaxMod led to a total lapse time of about 1.4s, while the equispaced approach led to 0.13s. The computation time of MaxMod is mostly due to the update of the covariance parameters (maximum likelihood estimation) each time a new knot or active variable is tried.

Thus, the execution time of MaxMod is larger than for the equispaced approach, for the same accuracy. However, the benefit of MaxMod is that it reduces the number of knots. This is crucial for the subsequent exploitation of the Gaussian process model. More precisely, the cost of numerical sampling of constrained conditional simulations of $f$, that enables to compute confidence intervals [38, 39, 42], increases with the number of knots. Here (in two dimensions) this cost is mild: obtaining 1000 constrained GP simulations required 0.08s with the knots provided by MaxMod, compared to 0.14s for the equispaced approach.<fnref n="2" /> However, the benefit of having less knots with MaxMod is key in higher dimension, as we illustrate below in dimension 5 (see Table 3). Finally, the equispaced approach becomes intractable when the dimension $D$ is larger than approximately 10, as the available numerical routines are unable to compute a mode with more than about a thousand ($2^{10}$) knots. This approach may be already inaccurate in the range 6 – 10 for $D$, as the common number of knots for the variables cannot exceed 3. In contrast, we will implement MaxMod up to dimension 20.

<footnote n="2">²In this paper, we considered the Hamiltonian Monte Carlo sampler proposed in [44] for simulating constrained GP realizations.</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 5: Evolution of the (normalized) bending energy $E_n$ for the example in Figure 3 with respect to the sample size $n = 10, \dots, 50$. The $E_n$ criterion is computed over a $10 \times 10$ equispaced grid.</figcaption></figure>

### 5.1.3 Impact of the sample size

In the context of expensive black box objective functions, we are interested in having accurate predictions with small sample sizes (number of observations of $f$). We assess the quality of the three modes for varying sample sizes $n = 10, \dots, 50$, considering again maximin LHDs. We use the optimal configurations of knots obtained in the experiment in Figure 4, that is $m = 15$ ($m_1 = 3, m_2 = 5$) for $\hat{Y}_{MaxMod,rect}$ and $\hat{Y}_{MaxMod}$ and $m = 16$ ($m_1 = m_2 = 4$) for $\hat{Y}_{square}$. Figure 5 shows that for small sample sizes, the three modes lead to small $E_n$ values, and that both $\hat{Y}_{MaxMod,rect}$ and $\hat{Y}_{MaxMod}$ outperform $\hat{Y}_{square}$. For this latter mode, we can see a local increase of the error from sample size 20 to 30, which may be explained by the randomness of the maximin LHDs. Hence, again, the prediction accuracy improves when considering the optimal design of knots of MaxMod.

## 5.2 Dimension reduction illustration

We now focus on the capability of MaxMod to perform dimension reduction. We apply the same stopping rule as in the previous test case, with a tolerance equal to $5 \times 10^{-3}$. We consider the target function:

$$f(x) = \sum_{i=1}^{d} \arctan \left( 5 \left[ 1 - \frac{i}{d+1} \right] x_i \right), \qquad (20)$$

with $x \in [0, 1]^d$. Note that $f$ is completely monotone exhibiting lesser growth rates as $i$ increases. In addition to $(x_1, \dots, x_d)$, we include $D - d$ virtual variables, indexed as $(x_{d+1}, \dots, x_D)$, which will compose the subset of inactive dimensions since $f$ does not depend on them. We consider $D \in \{5, 10, 15, 20\}$ and $d \in \{2, 3, 4, 5\}$. We also analyze the case where $d = D$. For each value of $D$, we evaluate $f$ at a maximin LHD with $n = 10 \times D$ points. For any possible combination of $D$ and $d$, we apply MaxMod expecting at properly finding the true $d$ active dimensions when $d < D$. When $d = D$, we expect that MaxMod will concentrate the computational budget on the most important input variables $x_1, \dots, x_{d^*}$ with $d^* < d$.

From Table 2, we observe that when $d < D$ MaxMod properly identifies the $d$ dimensions that are actually active, leading to small $E_n$ results when considering either $\hat{Y}_{MaxMod,rect}$ or $\hat{Y}_{MaxMod}$.<fnref n="3" /> When $d = D$, MaxMod successfully detects the most important variables $1, \dots, d^*$ with $d^* < D$ for $D \ge 10$. When the active dimension $d$ is larger or equal to 10, the function approximation problem is intrinsically more difficult. Thus MaxMod stops before being able to reach the stopping criterion, since the number of knots becomes too large for the numerical routines that compute the mode ($m > 1400$). In these cases, we fix the maximal number of iterations to 12, which nevertheless yields small values of $E_n$.

Overall, $\hat{Y}_{MaxMod}$ outperformed $\hat{Y}_{MaxMod,rect}$ due to its flexibility to freely allocate knots without being limited to equispaced designs. Finally, we remark that the mode $\hat{Y}_{square}$ is intractable here due to the large dimension $D$.

<footnote n="3">³In all the replicates in Table 2, MaxMod leads to small noise variance parameters $\hat{\tau}_{MaxMod}^2$ when $d < D$. The largest value, obtained for $D = 20$ and $d = 5$, is $\hat{\tau}_{MaxMod}^2 = 1.75 \times 10^{-3}$ ($\approx 0.4\%$ of the variance of the observations).</footnote>
<table><caption>Table 2: Performance of MaxMod for the example in Section 5.2. Results are shown for both $\hat{Y}_{\text{MaxMod,rect}}$ and $\hat{Y}_{\text{MaxMod}}$ considering $D \in \{5, 10, 15, 20\}$ and $d \in \{2, 3, 4, 5\}$. We also consider the case where $d = D$. The activated dimensions, number of one-dimensional knots per active dimension and $E_n$ results are displayed for any combination of $D$ and $d$. †The maximal number of iterations is fixed to 12 for the experiments where $d \ge 10$.</caption><thead><tr><th>D</th><th>d</th><th>active dim.</th><th>knots per dim.</th><th>E<sub>n</sub>(f, Y&#x0302;<sub>MaxMod,rect</sub>)</th><th>E<sub>n</sub>(f, Y&#x0302;<sub>MaxMod</sub>)</th></tr></thead><tbody><tr><td rowspan="4">5</td><td>2</td><td>(1, 2)</td><td>(7, 4)</td><td><b>2.59 &times; 10<sup>-5</sup></b></td><td>4.51 &times; 10<sup>-5</sup></td></tr><tr><td>3</td><td>(1, 2, 3)</td><td>(6, 6, 4)</td><td>6.42 &times; 10<sup>-4</sup></td><td><b>4.09 &times; 10<sup>-4</sup></b></td></tr><tr><td>4</td><td>(1, ..., 4)</td><td>(4, 4, 3, 2)</td><td><b>9.02 &times; 10<sup>-4</sup></b></td><td>9.05 &times; 10<sup>-4</sup></td></tr><tr><td>5</td><td>(1, ..., 5)</td><td>(3, 4, 4, 3, 2)</td><td><b>1.15 &times; 10<sup>-3</sup></b></td><td>1.19 &times; 10<sup>-3</sup></td></tr><tr><td rowspan="5">10</td><td>2</td><td>(1, 2)</td><td>(5, 3)</td><td><b>1.03 &times; 10<sup>-5</sup></b></td><td>2.78 &times; 10<sup>-5</sup></td></tr><tr><td>3</td><td>(1, 2, 3)</td><td>(5, 4, 3)</td><td>2.13 &times; 10<sup>-3</sup></td><td><b>1.79 &times; 10<sup>-3</sup></b></td></tr><tr><td>4</td><td>(1, ..., 4)</td><td>(5, 3, 3, 2)</td><td>3.10 &times; 10<sup>-4</sup></td><td><b>2.89 &times; 10<sup>-4</sup></b></td></tr><tr><td>5</td><td>(1, ..., 5)</td><td>(3, 4, 3, 3, 2)</td><td>7.40 &times; 10<sup>-4</sup></td><td><b>4.31 &times; 10<sup>-4</sup></b></td></tr><tr><td>10<sup>&dagger;</sup></td><td>(1, ..., 7)</td><td>(3, 3, 3, 2, 3, 3, 2)</td><td>4.44 &times; 10<sup>-3</sup></td><td><b>3.94 &times; 10<sup>-3</sup></b></td></tr><tr><td rowspan="5">15</td><td>2</td><td>(1, 2)</td><td>(4, 3)</td><td><b>1.73 &times; 10<sup>-4</sup></b></td><td>1.85 &times; 10<sup>-4</sup></td></tr><tr><td>3</td><td>(1, 2, 3)</td><td>(4, 3, 3)</td><td>2.13 &times; 10<sup>-4</sup></td><td><b>1.94 &times; 10<sup>-4</sup></b></td></tr><tr><td>4</td><td>(1, ..., 4)</td><td>(3, 3, 3, 2)</td><td>4.06 &times; 10<sup>-4</sup></td><td><b>1.94 &times; 10<sup>-4</sup></b></td></tr><tr><td>5</td><td>(1, ..., 5)</td><td>(3, ..., 3, 2)</td><td>2.26 &times; 10<sup>-4</sup></td><td><b>9.29 &times; 10<sup>-5</sup></b></td></tr><tr><td>15<sup>&dagger;</sup></td><td>(1, ..., 6)</td><td>(3, ..., 3)</td><td>6.22 &times; 10<sup>-3</sup></td><td><b>5.93 &times; 10<sup>-3</sup></b></td></tr><tr><td rowspan="5">20</td><td>2</td><td>(1, 2)</td><td>(5, 3)</td><td>9.88 &times; 10<sup>-5</sup></td><td><b>9.37 &times; 10<sup>-5</sup></b></td></tr><tr><td>3</td><td>(1, 2, 3)</td><td>(4, 4, 3)</td><td><b>1.40 &times; 10<sup>-4</sup></b></td><td><b>1.40 &times; 10<sup>-4</sup></b></td></tr><tr><td>4</td><td>(1, ..., 4)</td><td>(4, 3, 3, 3)</td><td>3.48 &times; 10<sup>-4</sup></td><td><b>1.97 &times; 10<sup>-4</sup></b></td></tr><tr><td>5</td><td>(1, ..., 5)</td><td>(3, ..., 3, 2)</td><td>5.60 &times; 10<sup>-4</sup></td><td><b>2.83 &times; 10<sup>-4</sup></b></td></tr><tr><td>20<sup>&dagger;</sup></td><td>(1, ..., 7)</td><td>(2, 3, ..., 3, 2)</td><td>5.58 &times; 10<sup>-3</sup></td><td><b>5.55 &times; 10<sup>-3</sup></b></td></tr></tbody></table>

## 5.3 Coastal flooding application in 5D

We consider the 5D coastal flood application studied in [2, 38], available in the R package *profExtrema* [1]. In the past, a flood event at the *Boucholeurs* area (La Rochelle, France) was induced by an overflow on the Atlantic ocean caused by the Xynthia storm in 2010. To prevent adverse coastal flood events, such as the one led by the Xynthia storm, accurate forecast and early-warning systems (see, e.g., [2, 38, 48] for GP-based ones) are required.

The dataset contains 200 observations of the flooded area ($A_{\text{flood}}[m^2]$) driven by five offshore forcing conditions (inputs) at the *Boucholeurs* area: tide ($T[m]$), surge ($S[m]$), the phase difference ($\phi$, hours) between the surge peak and the high tide, the time duration of the raising part ($t_-$, hours) and the falling part ($t_+$, hours) of the (triangular) surge signal. In particular, it is known that $A_{\text{flood}}$ increases as $T$ and $S$ increase. According to [2, 38], while the contribution of $T$, $S$, $t_-$ and $t_+$ are almost linear, $A_{\text{flood}}$ exhibits a higher variation across $\phi$. Thus, in our experiments, we may expect that MaxMod properly concentrate one-dimensional knots across $\phi$ rather than across the other dimensions. As suggested in [38], we consider here $Y := \log_{10}(A_{\text{flood}})$ as the output variable.

As shown in [38], enforcing a GP-based coastal emulator to both positivity and monotonicity (with respect to $T$ and $S$) constraints leads to a more reliable prediction. There, the number of one-dimensional knots per dimension has been manually fixed looking for a trade-off between the computational cost and the quality of resolution of the constrained GP. Here, we aim at applying MaxMod and comparing the $E_n$ results to those from the knots in [38]. As in Sections 5.1 and 5.2, we compute the modes $\hat{Y}_{\text{square}}$, $\hat{Y}_{\text{MaxMod,rect}}$ and $Y_{\text{MaxMod}}$. In addition, we compute the mode resulting from the knots in [38], denoted as $\hat{Y}_*$. Since the target function $Y$ is actually unknown, the bending energy $E_n$ is computed over the available 200 observations: $E_n(Y, \hat{Y}) = \sum_{i=1}^{200} (Y_i - \hat{Y}_i)^2 / \sum_{i=1}^{200} Y_i^2$. Figure 6 shows that MaxMod results in a total of $m = 4 \times 3 \times 6 \times 3 \times 2 = 432$ knots, leading to $E_n(Y, \hat{Y}_{\text{MaxMod,rect}}) = 9.17 \times 10^{-3}$ and $E_n(Y, \hat{Y}_{\text{MaxMod}}) = 8.81 \times 10^{-3}$. These results are comparable to those from $\hat{Y}_*$, $E_n(Y, \hat{Y}_*) = 8.72 \times 10^{-3}$ with $m = 4 \times 4 \times 5 \times 3 \times 3 = 720$, and from $\hat{Y}_{\text{square}}$, $E_n(Y, \hat{Y}_{\text{square}}) = 8.72 \times 10^{-3}$ with $m = 4^5 = 1024$. Moreover, as expected, MaxMod concentrated the computational budget on the input $\phi$ (6 one-dimensional knots) rather than the other ones
<figure><img src="image_7.png" /><figcaption>Figure 6: Evolution of the bending energy $E_n$ criterion for the coastal flooding application in Section 5.3. The panel description is the same as in Figure 4. The $E_n$ value associated to the mode $\hat{Y}_*$, yielded by using the configuration of the knots suggested in [38], is displayed by a purple asterisks at $m = 720$.</figcaption></figure>

<table><caption>Table 3: CPU lapse times for the coastal flooding application in Section 5.3.</caption><thead><tr><th rowspan="2">Approach</th><th rowspan="2">m</th><th rowspan="2">E<sub>n</sub>(Y, Ŷ<sub>Y</sub>)<br>[1 × 10<sup>-3</sup>]</th><th colspan="3">CPU time [s]</th></tr><tr><th>Training step</th><th>Computation<br>of Ŷ<sub>Y</sub></th><th>Sampling step<br>with 100 realizations</th></tr></thead><tbody><tr><td>ŷ<sub>square</sub></td><td>1024</td><td>8.72</td><td>49.1</td><td>8.03</td><td>non converged after 1 day</td></tr><tr><td>ŷ<sub>MaxMod</sub></td><td>432</td><td>8.81</td><td>949.5</td><td>0.58</td><td>108.72</td></tr></tbody></table>

(4,3,3,2) since Y varies the most across $\phi$. In Figure 6, it appears that $\hat{Y}_{MaxMod}$ is the most efficient, needing a smaller number of knots than $\hat{Y}_{square}$, $\hat{Y}_{MaxMod,rect}$ and $\hat{Y}_*$ to reach a given value of $E_n$. In terms of noise variance, MaxMod estimates $\hat{\sigma}_{MaxMod}^2 = 6.16 \times 10^{-2}$, a small (but non negligible) value equivalent to $\approx 7.6\%$ of the variance of the observations. Our interpretation for this small value is that it accounts for possible numerical instabilities of the computer code. Furthermore it improves the accuracy of the mode function and speeds up its computation by making the inequality constraints easier to satisfy (see [38] for further discussions).

Table 3 shows the CPU times required for training the GP models from $\hat{Y}_{square}$ (covariance parameter estimation) and $\hat{Y}_{MaxMod}$ (Algorithm 2 including covariance parameter estimation). It also shows the elapsed times related to the computation of the modes $\hat{Y}_{square}$ and $\hat{Y}_{MaxMod}$ (with given covariance parameters), and to the sampling of 100 conditional realizations of the resulting constrained processes. The training time is smaller for $\hat{Y}_{square}$. The computation of $\hat{Y}$ is faster for MaxMod, but for both MaxMod and the equispaced approach, the computation time for $\hat{Y}$ is small compared to the training time. With the 432 knots from MaxMod, the 100 conditional realizations are obtained in less than 2 minutes. In contrast, with the 1024 knots from $\hat{Y}_{square}$, the Hamiltonian Monte Carlo routine for sampling the conditional distribution proposed in [44] does not converge, because of the overly high dimension. This highlights the benefit of MaxMod in real applications: it requires an acceptable increase of the (offline) training time in order to reduce the constrained GP model dimension, which enhances its subsequent exploitation (conditional sampling).

# 6 Conclusion

This paper introduces MaxMod, that sequentially inserts one-dimensional knots or adds active variables to a constrained GP regression model. This algorithm results in the first constrained GP model that at the same time satisfies the constraints everywhere and that is not restricted to small dimensional cases in practice.

A proof of convergence, for a fixed dataset and as the number of iterations goes to infinity, guarantees that, despite its sequential nature, MaxMod globally converges to an optimal infinite dimensional model. In establishing this convergence, the notion of a multi-affine extension is
constructed. Furthermore, the convergence of finite-dimensional GP models is shown in settings where the multi-dimensional knots are not dense in the input domain, thereby extending the recent literature. This construction, and this latter convergence result, may be of special and independent interest, together with the corresponding proof techniques.

In Section 5, with simulated and real data, it is confirmed numerically that MaxMod is tractable and accurate (at least in dimension $D = 20$), and typically needs less multi-dimensional knots than the other state-of-the-art constrained GP models. We demonstrate the strong benefit of having fewer knots when subsequently computing confidence intervals from constrained GP models. These numerical examples also indicate that MaxMod successfully detects the active variables when the effective dimension $d$ is small. Even when $d = D$, MaxMod was able to detect the most important variables.

## Acknowledgments

We are grateful to N. Durrande for constructive discussions. This research was supported by the Chair in Applied Mathematics OQUAIDO<fnref n="4" />, the ANR<fnref n="5" />, and the Isaac Newton Institute for Mathematical Sciences, Cambridge.<fnref n="6" /> Part of this research was carried out when the second author was affiliated to IMT and BRGM. We are indebted to the referees, for their constructive suggestions that greatly improved the manuscript.

## A Expression of the linear inequality constraints in Section 2.3.2 for boundedness, monotonicity and componentwise convexity

Let us provide expressions of $M(C_J)$ and of $v(C_J)$ such that (7) holds when $C_D$ is given by (3), (4) and (5). When $C_D$ is given by (3), we define $M(C_J)$ and $v(C_J)$ as follows. The matrix $M(C_J)$ has $2|C_S|$ rows and thus $v(C_J)$ is of size $2|C_S|$. The first $|C_S|$ rows and values of $v(C_J)$ are indexed by $\underline{\ell} \in C_S$ and defined by, for $\ell' \in C_S$,

$$M(C_J)_{\underline{\ell}, \ell'} = -\mathbf{1}_{\underline{\ell} = \ell'} \quad \text{and} \quad v(C_J)_{\underline{\ell}} = -a. \qquad (21)$$

The last $|C_S|$ rows and values of $v(C_J)$ are indexed by $\underline{\ell} \in C_S$ and defined by, for $\ell' \in C_S$,

$$M(C_J)|_{C_S|+\underline{\ell}, \ell'} = \mathbf{1}_{\underline{\ell} = \ell'} \quad \text{and} \quad v(C_J)|_{C_S|+\underline{\ell}} = b, \qquad (22)$$

using the slight abuse of notation $|C_S| + \underline{\ell}$ to denote the last $|C_S|$ rows and values. The fact that (7) holds with (21) and (22) can be simply shown to follow from the proof of Lemma 3 and from (19). Similar equations were also stated, in the previous references [39, 42].

When $C$ is given by (4), we define $M(C_J)$ and $v(C_J)$ as follows. The matrix $M(C_J)$ is composed of $d = |\mathcal{J}|$ vertically stacked matrices $M^{(1)}(C_J), \dots, M^{(d)}(C_J)$ and $v(C_J)$ is composed of $d$ vertically stacked column vectors $v^{(1)}(C_J), \dots, v^{(d)}(C_J)$. For $i = 1, \dots, d$, the matrix $M^{(i)}(C_J)$ has $(m_{S_{a_i}} - 1) \prod_{j=1, \dots, d, j \neq i} m_{S_{a_j}}$ rows and $v^{(i)}(C_J)$ has the same number of components. The rows and values of $v^{(i)}(C_J)$ are indexed by the multi-indices $\underline{\ell} = (\ell_{a_1}, \dots, \ell_{a_d}) \in (\prod_{j=1}^{i-1} \{1, \dots, m_{S_{a_j}}\}) \times \{2, \dots, m_{S_{a_i}}\} \times (\prod_{j=i+1}^{d} \{1, \dots, m_{S_{a_j}}\})$ and defined by

$$M^{(i)}(C_J)_{(\ell_{a_1}, \dots, \ell_{a_d}), (\ell'_{a_1}, \dots, \ell'_{a_d})} = -\mathbf{1}_{\ell_{a_j} = \ell'_{a_j}} \quad \text{for } j \neq i, \ell_{a_i} = \ell'_{a_i} + \mathbf{1}_{\ell_{a_j} = \ell'_{a_j}} \quad \text{for } j \neq i, \ell_{a_i} - 1 = \ell'_{a_i} \qquad (23)$$
and $v^{(i)}(C_J)_{(\ell_{a_1}, \dots, \ell_{a_d})} = 0$.

Equation (23) means that the set of values on the tensorized grid of $d$-dimensional knots is component-wise non-decreasing. It can be shown from the proof of Lemma 3 and from (19) that $M(C_J)$ and $v(C_J)$ given by (23) imply (7). This was noticed in the previous references [39, 42].

Finally, when $C$ is given by (5), we define $M(C_J)$ and $v(C_J)$ as follows. The matrix $M(C_J)$ is composed of $d$ vertically stacked matrices $M^{(1)}(C_J), \dots, M^{(d)}(C_J)$ and $v(C_J)$ is composed of $d$

<footnote n="4">⁴OQUAIDO gathers partners in technological research (BRGM, CEA, IFPEN, IRSN, Safran, Storengy) and academia (CNRS, Ecole Centrale de Lyon, Mines Saint-Etienne, Univ. of Grenoble, Univ. of Nice, Univ. of Toulouse) around advanced methods for Computer Experiments.</footnote>

<footnote n="5">⁵French National Research Agency, under the RISCOPE project (ANR-16-CE04-0011).</footnote>

<footnote n="6">⁶UNQ program (EPSRC grantEP/K032208/1).</footnote>
vertically stacked column vectors $v^{(1)}(\mathcal{C}_J),...,v^{(d)}(\mathcal{C}_J)$. For $i=1,...,d$, the matrix $M^{(i)}(\mathcal{C}_J)$ has $(m_{S_{a_i}} - 2) \prod_{j=1,...,d,j \neq i} m_{S_{a_j}}$ rows and $v^{(i)}(\mathcal{C}_J)$ has the same number of components. The rows and values of $v^{(i)}(\mathcal{C}_J)$ are indexed by the multi-indices $\ell = (\ell_{a_1},...,\ell_{a_d}) \in (\prod_{j=1}^{i-1} \{1,...,m_{S_{a_j}}\}) \times \{3,...,m_{S_{a_i}}\} \times (\prod_{j=i+1}^{d} \{1,...,m_{S_{a_j}}\})$ and defined by

$$ M^{(i)}(\mathcal{C}_J)_{(\ell_{a_1},...,\ell_{a_d}), (\ell'_{a_1},...,\ell'_{a_d})} = -\mathbf{1}_{\ell_{a_j}=\ell'_{a_j}} \quad \text{for } j \neq i, \ell_{a_i}=\ell'_{a_i} + 2\mathbf{1}_{\ell_{a_j}=\ell'_{a_j}} \quad \text{for } j \neq i, \ell_{a_i}-1=\ell'_{a_i} - \mathbf{1}_{\ell_{a_j}=\ell'_{a_j}} \quad \text{for } j \neq i, \ell_{a_i}-2=\ell'_{a_i} \quad (24) $$

and $v^{(i)}(\mathcal{C}_J)_{(\ell_{a_1},...,\ell_{a_d})} = 0$.

Equation (24) means that the set of values on the tensorized grid of $d$-dimensional knots is component-wise convex. It follows from the proof of Lemma 3 and from (19) that $M(\mathcal{C}_J)$ and $v(\mathcal{C}_J)$ given by (24) imply (7). To the best of our knowledge, this was only shown in dimension one, in earlier references.

**Remark 3.** The equation (24) is defined only for $i \in \{1,...,d\}$ such that $m_{S_{a_i}} \ge 3$. If $i$ is such that $m_{S_{a_i}} = 2$, then the matrix $M^{(i)}(\mathcal{C}_J)$ and the vector $v^{(i)}(\mathcal{C}_J)$ can be removed from $M(\mathcal{C}_J)$ and $v(\mathcal{C}_J)$. Indeed, when $m_{S_{a_i}} = 2$, the function $\pi_S(\xi_J)$ is linear (thus convex) along the dimension $a_i$.

# B Computing the $L^2$ difference between modes in practice

At step $m$ of the MaxMod algorithm, the current mode function is $\hat{Y}_{\mathcal{J}_m, S^{(m)}} : [0,1]^{d_m} \to \mathbb{R}$, with $d_m = |\mathcal{J}_m|$ and with

$$ \hat{Y}_{\mathcal{J}_m, S^{(m)}} = \sum_{\ell \in \mathcal{L}_{S^{(m)}}} (\hat{\alpha}_{\mathcal{J}_m, S^{(m)}})_{\ell} \phi_{\ell}^{(S^{(m)})}. \quad (25) $$

To simplify the exposition in Section B, we let $\mathcal{J}_m = \mathcal{J}$, and we consider the case where $\mathcal{J} = \{1,...,d\}$. The results of Section B can then be immediately extended to a general set $\mathcal{J}$, by permuting indices. We also let $S^{(m)} = S = (S_1,...,S_d)$.

## B.1 $L^2$ difference for a knot insertion to an active variable

Consider the active variable 1, and a new knot $t \in [0,1]\setminus S_1^{(m)}$ (for the variable 1). Again, the choice of the active variable 1 enables to simplify the exposition and the results of Section B.1 can be immediately extended to a general active variable $i \in \{1,...,d\}$.

The candidate mode function is

$$ \hat{Y}_{\mathcal{J}, S \cup_1 t} = \sum_{\ell \in \mathcal{L}_{S \cup_1 t}} (\hat{\alpha}_{\mathcal{J}, S \cup_1 t})_{\ell} \phi_{\ell}^{(S \cup_1 t)}. \quad (26) $$

For $j = 1,...,d$, write $m_j = m_{S_j}$. With this notation, $S_j = \{t_{(0)}^{(S_j)},...,t_{(m_j+1)}^{(S_j)}\}$ with the ordered knots $t_{(0)}^{(S_j)} < \dots < t_{(m_j+1)}^{(S_j)}$.

Let $S' = S \cup_1 t = (S'_1, ..., S'_d)$. We have $S'_j = S_j$ for $j \in \{2,...,d\}$. Let $\nu \in \{1,...,m_1-1\}$ be such that $t_{(\nu)}^{S_1} < t < t_{(\nu+1)}^{S_1}$. Then we have $S'_1 = \{t_{(0)}^{(S'_1)},...,t_{(m_1+2)}^{(S'_1)}\}$ with the ordered knots $t_{(0)}^{(S'_1)} < \dots < t_{(m_1+2)}^{(S'_1)}$ with

$$ (t_{(0)}^{(S'_1)},...,t_{(m_1+2)}^{(S'_1)}) = (t_{(0)}^{(S_1)},...,t_{(\nu)}^{(S_1)},t_{(\nu+1)}^{(S_1)},...,t_{(m_1+1)}^{(S_1)}). $$

Then the next proposition provides a computationally efficient formula for the $L^2$ difference between modes. To understand Proposition 3, note that the difference between the current mode function (25) and the new mode function (26) can be expressed as a weighted sum of the basis functions of the new (refined) subdivision $S'$. In Proposition 3, (27), (28) and (29) first provide the expressions of the corresponding coefficients. Then, it just suffices to compute the matrix of $L^2$ inner products between the $d$-dimensional basis functions of $S'$. This matrix is given in (30). It is the tensor product of the $d$ matrices of $L^2$ inner products between the one-dimensional basis functions of the $d$ one-dimensional subdivisions $S'_1, ..., S'_d$. These latter matrices are provided in (31).
**Proposition 3** ($L^2$ difference for knot insertion to an active variable). Define for $\underline{\ell} = (\ell_1, \dots, \ell_d) \in \mathcal{L}_{S'}$,

$$ \beta_{\underline{\ell}} = (\hat{\alpha}_{\mathcal{J},S})_{\underline{\ell}} - (\hat{\alpha}_{\mathcal{J},S'})_{\underline{\ell}} \quad (27) $$

if $\ell_1 \in \{1, \dots, \nu\}$,

$$ \beta_{\underline{\ell}} = (\hat{\alpha}_{\mathcal{J},S})(\underline{\ell}, \ell_2, \dots, \ell_d) \frac{t_{(\nu+1)}^{(S_1)} - t}{t_{(\nu+1)}^{(S_1)} - t_{(\nu)}^{(S_1)}} + (\hat{\alpha}_{\mathcal{J},S})(\underline{\ell}+1, \ell_2, \dots, \ell_d) \frac{t - t_{(\nu)}^{(S_1)}}{t_{(\nu+1)}^{(S_1)} - t_{(\nu)}^{(S_1)}} - (\hat{\alpha}_{\mathcal{J},S'})_{\underline{\ell}}, \quad (28) $$

if $\ell_1 = \nu + 1$ and

$$ \beta_{\underline{\ell}} = (\hat{\alpha}_{\mathcal{J},S})(\underline{\ell}-1, \ell_2, \dots, \ell_d) - (\hat{\alpha}_{\mathcal{J},S'})_{\underline{\ell}} \quad (29) $$

if $\ell_1 \in \{\nu + 2, \dots, m_1 + 1\}$.

Define for $\underline{\ell} = (\ell_1, \dots, \ell_d) \in \mathcal{L}_S$, $\underline{\ell}' = (\ell'_1, \dots, \ell'_d) \in \mathcal{L}_{S'}$,

$$ \Psi_{\underline{\ell}, \underline{\ell}'} = \begin{cases} \prod_{j=1, \dots, d} \psi_{\ell_j, \ell_j'}^{(j)} & \text{if } |\ell_j - \ell_j'| \le 1 \text{ for } j = 1, \dots, d, \\ 0 & \text{else} \end{cases}, \quad (30) $$

with, when $|\ell_j - \ell_j'| \le 1$,

$$ \psi_{\ell_j, \ell_j'}^{(j)} = \begin{cases} \frac{t_{(\ell_j+1)}^{(S')} - t_{(\ell_j)}^{(S')}}{3} & \text{if } \ell_j = \ell_j' = 1 \\ \frac{t_{(\ell_j)}^{(S')} - t_{(\ell_j-1)}^{(S')}}{3} & \text{if } \ell_j = \ell_j' = m_j + 1 \\ \frac{t_{(\ell_j+1)}^{(S')} - t_{(\ell_j-1)}^{(S')}}{3} & \text{if } \ell_j = \ell_j' \in \{2, \dots, m_j\} \\ \frac{t_{(\ell_j+1)}^{(S')} - t_{(\ell_j)}^{(S')}}{6} & \text{if } \ell_j' = \ell_j + 1 \\ \frac{t_{(\ell_j)}^{(S')} - t_{(\ell_j-1)}^{(S')}}{6} & \text{if } \ell_j' = \ell_j - 1 \end{cases}. \quad (31) $$

Then, we have

$$ \int_{[0,1]^d} (\hat{Y}_{\mathcal{J},S}(x) - \hat{Y}_{\mathcal{J},S \cup t}(x))^2 dx = \sum_{\underline{\ell} \in \mathcal{L}_{S'}, \underline{\ell}' \in \mathcal{L}_{S'}} \beta_{\underline{\ell}} \Psi_{\underline{\ell}, \underline{\ell}'} \beta_{\underline{\ell}'} \quad (32) $$

Consider that the coefficients in (25) and (26) have been computed by optimization of a quadratic function with linear inequality constraints, see Section 2.3. Then, from the above proposition, the $L^2$ difference between modes can be obtained by the explicit quadratic form (32). Furthermore, the matrix $\Psi$ defining this quadratic form is the tensor product of $d$ banded matrices $\psi^{(1)}, \dots, \psi^{(d)}$. Hence, the computational cost is eventually linear in the number of multidimensional knots.

## B.2 $L^2$ difference when a new active variable is added

Consider the new variable $d+1$. Again, the choice of the new variable $d+1$ enables to simplify the exposition and the results of Section B.2 can be immediately extended to a general new variable $i \in \{d+1, \dots, D\}$. The candidate mode function is

$$ \hat{Y}_{\mathcal{J} \cup \{d+1\}, S+(d+1)} = \sum_{\underline{\ell} \in \mathcal{L}_{S+(d+1)}} (\hat{\alpha}_{\mathcal{J} \cup \{d+1\}, S+(d+1)})_{\underline{\ell}} \phi_{\underline{\ell}}^{(S+(d+1))}. \quad (33) $$

Then the next proposition provides a computationally efficient formula for the $L^2$ difference between modes.

**Proposition 4** ($L^2$ difference when a new active variable is added). Let $\underline{\ell} \in \mathcal{L}_{S+(d+1)}$ and write $\underline{\ell}$ of the form $(\underline{\ell}, \ell_{d+1})$ with $\underline{\ell} \in \mathcal{L}_S$ and $\ell_{d+1} \in \{1, 2\}$. Then let

$$ \beta_{\underline{\ell}} = (\hat{\alpha}_{\mathcal{J},S})_{\underline{\ell}} - (\hat{\alpha}_{\mathcal{J} \cup \{d+1\}, S+(d+1)})_{\underline{\ell}}. $$

For $\underline{\ell}, \underline{\ell}' \in \mathcal{L}_{S+(d+1)}$, define $\Psi_{\underline{\ell}, \underline{\ell}'}$ as in Proposition 3, but with $S'$ there replaced by $S + (d+1)$. Then we have

$$ \int_{[0,1]^{d+1}} (\hat{Y}_{\mathcal{J},S}(x) - \hat{Y}_{\mathcal{J} \cup \{d+1\}, S+(d+1)}(x))^2 dx = \sum_{\underline{\ell} \in \mathcal{L}_{S+(d+1)}, \underline{\ell}' \in \mathcal{L}_{S+(d+1)}} \beta_{\underline{\ell}} \Psi_{\underline{\ell}, \underline{\ell}'} \beta_{\underline{\ell}'} \quad (34) $$

The discussion of Proposition 4 is the same as for Proposition 3.
C Technical conditions for Theorem 3

As for Theorem 2, for each $\mathcal{J} \subseteq \{1, \dots, D\}$, we consider a set of functions in $C([0, 1]^{|\mathcal{J}|}, \mathbb{R})$ satisfying strict inequalities: $\dot{\mathcal{C}}_{\mathcal{J}} \subseteq \mathcal{C}_{\mathcal{J}}$. We recall that the MaxMod algorithm is initialized with $\mathcal{J}_0 \subseteq \{1, \dots, D\}$ and $S^{(0)} \in \mathcal{S}$.

Conditions 2, 3 and 4 are replaced by the following conditions.

**Condition 7** (extension of Condition 2). Consider $\mathcal{J} \supseteq \mathcal{J}_0$ and $S \in \mathcal{S}_{\mathcal{J}}$ such that, for $j \in \mathcal{J}_0$, we have $S_j \supseteq S_j^{(0)}$. Then the set

$$ \{\alpha \in A_S; Y_{S,\alpha} \in \mathcal{I}_{X_{\mathcal{J},y^{(n)}}} \cap \dot{\mathcal{C}}_{\mathcal{J}}\} $$

is non-empty.

Condition 7 means that whenever the set of active variables has been increased compared to $\mathcal{J}_0$ and knots have been inserted to the variables in $\mathcal{J}_0$ or to the new active variables, then it is possible to find an interpolating function that satisfies the strict inequalities.

**Condition 8** (extension of Condition 3). We have

$$ \left\{ \left( Y_{S^{(0)},\alpha} \left( x_{\mathcal{J}_0}^{(i)} \right) \right)_{i=1,\dots,n} ; \alpha \in S^{(0)} \right\} = \mathbb{R}^n. $$

For $\mathcal{J} \subseteq \{1, \dots, D\}$, we let $\mathcal{H}_{\mathcal{J}}$ be the RKHS of $k_{\mathcal{J}}$ with norm $||\cdot||_{\mathcal{H}_{\mathcal{J}}}$ (recall that $\mathcal{H}_{\mathcal{J}}$ is a set of functions from $[0, 1]^{|\mathcal{J}|}$ to $\mathbb{R}$).

**Condition 9** (extension of Condition 4). For all $\mathcal{J} \subseteq \{1, \dots, D\}$, for all $h \in \dot{\mathcal{C}}_{\mathcal{J}}$, for all $r \in \mathbb{N}$, $U = (u_1, \dots, u_r) \in ([0, 1]^{|\mathcal{J}|})^r$, with $u_1, \dots, u_r$ two-by-two distinct, letting $v^{(r)} = (h(u_1), \dots, h(u_r)) \in \mathbb{R}^r$, the set $\text{int}_{||\cdot||_{\mathcal{H}_{\mathcal{J}}}}(\mathcal{H}_{\mathcal{J}} \cap \mathcal{C}_{\mathcal{J}}) \cap \mathcal{I}_{U,v^{(r)}}$ is non-empty.

Consider the case where $\mathcal{C}_D$ is given by (3), (4) and (5) (boundedness, monotonicity and convexity). Then the same discussions of Conditions 2 and 4 of Theorem 2 apply to Conditions 7 and 9 here. In particular, when $\mathcal{C}_D$ is given by (3), we choose $\dot{\mathcal{C}}_{\mathcal{J}}$ as the set of functions in $C([0, 1]^{|\mathcal{J}|}, \mathbb{R})$ that are strictly between $a$ and $b$. When $\mathcal{C}_D$ is given by (4), we choose $\dot{\mathcal{C}}_{\mathcal{J}}$ as the set of functions in $C([0, 1]^{|\mathcal{J}|}, \mathbb{R})$ that are strictly increasing. When $\mathcal{C}_D$ is given by (5), we choose $\dot{\mathcal{C}}_{\mathcal{J}}$ as the set of functions in $C([0, 1]^{|\mathcal{J}|}, \mathbb{R})$ which one-dimensional cuts are strictly convex. With these choices, Conditions 7 and 9 are mild and mean that there are sufficiently many knots and sufficiently many active variables at the beginning of the MaxMod algorithm and that the Hilbert space $\mathcal{H}_{\mathcal{J}}$ is sufficiently rich for $\mathcal{J} \supseteq \mathcal{J}_0$. Finally, Lemma 2 can be straightforwardly adapted.

Condition 5 is replaced by the following one.

**Condition 10** (extension of Condition 5). For all $\mathcal{J} \subseteq \{1, \dots, D\}$, for all subdivision $S \in \mathcal{S}_{\mathcal{J}}$, $\pi_S(\mathcal{C}_{\mathcal{J}}) \subseteq \mathcal{C}_{\mathcal{J}}$.

As for Theorem 2, Condition 10 can be shown to hold for boundedness, monotonicity and input-wise convexity constraints. Finally, Condition 6 is replaced by the following one.

**Condition 11** (extension of Condition 6). For any $\mathcal{J} = (a_1, \dots, a_d) \subseteq \{1, \dots, D\}$, for a closed set $F = F_{a_1} \times \dots \times F_{a_d} \subseteq [0, 1]^d$ with $0, 1 \in F_{a_j}$ for $j = 1, \dots, d$, let us define $P_{F \to [0,1]^d}$ as in Proposition 1. Then, for $f \in \mathcal{C}_{\mathcal{J}}$, with $f|_F$ the restriction of $f$ to $F$, we have $P_{F \to [0,1]^d}(f|_F) \in \mathcal{C}_{\mathcal{J}}$.

Remark that Lemma 3 can be extended straightforwardly.

D Proofs

**Proof of Lemma 1.** Since the reward $R_{\mathcal{J}_m, S^{(m)}}(i, t)$ is smaller than $\max(\Delta, \Delta')$, it is sufficient to show that, for $i \in \{1, \dots, D\}$,

$$ \begin{array}{ll} \displaystyle \sup_{t \in [0,1],} & I_{\mathcal{J}_m, S^{(m)}}(i,t) < \infty. \tag{35} \\ \text{s.t.} & d(t, S_m^{(m)}) \ge b_m \\ & \text{if } i \in \mathcal{J}_m \end{array} $$

If $i \notin \mathcal{J}_m$, $I_{\mathcal{J}_m, S^{(m)}}(i, t)$ in (35) does not depend on $t$ (since the knots for the variable $i$, after this variable has been added, are $\{-1, 0, 1, 2\}$ independently on $t$). Thus the sup in (35) is finite.
Consider now $i \in \mathcal{J}_m$. Let $E_{i,m}$ be the set of $t \in [0,1]$ such that $d(t, S_i^{(m)}) \ge b_m$. By continuity of $d(\cdot, S_i^{(m)})$, notice that $E_{i,m}$ is compact.

Consider the function $\hat{Y}_{\mathcal{J}_m, S^{(m)}}$. We can write it in the thinner space obtained by inserting the knot $t$ at coordinate $i$, in the form $Y_{S^{(m)} \cup_i t, \alpha_{m,t}}$ for some $\alpha_{m,t} \in A_{S^{(m)} \cup_i t}$.

The coefficients $\alpha_{m,t}$ are equal to the values of the function $\hat{Y}_{\mathcal{J}_m, S^{(m)}}$ at the (multidimensional) knots in $S^{(m)} \cup_i t$, as remarked just after (1). Since $\hat{Y}_{\mathcal{J}_m, S^{(m)}}$ does not depend on $t$ and is continuous on $[0,1]^{|\mathcal{J}_m|}$, it follows that $\alpha_{m,t}$ is bounded with respect to $t \in E_{i,m}$.

Furthermore, let $\lambda_{\max}(M)$ and $\lambda_{\min}(M)$ be the largest and smallest eigenvalues of a matrix $M$. By definition of $\hat{\alpha}_{\mathcal{J}_m, S^{(m)} \cup_i t}$, and since $\hat{Y}_{\mathcal{J}_m, S^{(m)}}$ belongs to $\mathcal{I}_{X(\mathcal{J}_m), y^{(n)}} \cap \mathcal{C}_{\mathcal{J}_m}$, we have:

$$ \begin{aligned} & \hat{\alpha}_{\mathcal{J}_m, S^{(m)} \cup_i t}^{\top} k_{\mathcal{J}_m} (S^{(m)} \cup_i t, S^{(m)} \cup_i t)^{-1} \hat{\alpha}_{\mathcal{J}_m, S^{(m)} \cup_i t} \\ & \le \alpha_{m,t}^{\top} k_{\mathcal{J}_m} (S^{(m)} \cup_i t, S^{(m)} \cup_i t)^{-1} \alpha_{m,t} \\ & \le \frac{\|\alpha_{m,t}\|^2}{\lambda_{\min} (k_{\mathcal{J}_m} (S^{(m)} \cup_i t, S^{(m)} \cup_i t))}. \end{aligned} $$

Furthermore,

$$ \begin{aligned} & \hat{\alpha}_{\mathcal{J}_m, S^{(m)} \cup_i t}^{\top} k_{\mathcal{J}_m} (S^{(m)} \cup_i t, S^{(m)} \cup_i t)^{-1} \hat{\alpha}_{\mathcal{J}_m, S^{(m)} \cup_i t} \\ & \ge \frac{\|\hat{\alpha}_{\mathcal{J}_m, S^{(m)} \cup_i t}\|^2}{\lambda_{\max} (k_{\mathcal{J}_m} (S^{(m)} \cup_i t, S^{(m)} \cup_i t))}. \end{aligned} $$

Hence, we have that

$$ \sup_{t \in E_{i,m}} \|\hat{\alpha}_{\mathcal{J}_m, S^{(m)} \cup_i t}\|^2 \le \sup_{t \in E_{i,m}} \frac{\lambda_{\max} (k_{\mathcal{J}_m} (S^{(m)} \cup_i t, S^{(m)} \cup_i t))}{\lambda_{\min} (k_{\mathcal{J}_m} (S^{(m)} \cup_i t, S^{(m)} \cup_i t))} \sup_{t \in E_{i,m}} \|\alpha_{m,t}\|^2. $$

By assumption, $k_{\mathcal{J}_m}$ is continuous and the matrix $k_{\mathcal{J}_m}(S^{(m)} \cup_i t, S^{(m)} \cup_i t)$ is invertible for all $t$ in the compact set $E_{i,m}$. Thus, the ratio of eigenvalues above is bounded. Hence $\|\hat{\alpha}_{\mathcal{J}_m, S^{(m)} \cup_i t}\|$ is bounded with respect to $t \in E_{i,m}$ and thus also the supremum of the function $\hat{Y}_{\mathcal{J}_m, S^{(m)} \cup_i t}$ is bounded with respect to $t \in E_{i,m}$. Hence the sup in (35) is indeed finite, which concludes the proof. $\square$

**Proof of Proposition 3.** The current mode function

$$ \hat{Y}_{\mathcal{J}, S} = \sum_{\underline{\ell} \in \mathcal{L}_S} (\hat{\alpha}_{\mathcal{J}, S})_{\underline{\ell}} \phi_{\underline{\ell}}^{(S)} \quad (36) $$

can be expressed of the form

$$ \hat{Y}_{\mathcal{J}, S} = \sum_{\underline{\ell} \in \mathcal{L}_{S'}} \gamma_{\underline{\ell}} \phi_{\underline{\ell}}^{(S')} $$

Let us express $\gamma_{\underline{\ell}}$. Let $\underline{\ell} = (\ell_1, \dots, \ell_d) \in \mathcal{L}_{S'}$. First, if $\ell_1 \in \{1, \dots, \nu\}$, we have

$$ \gamma_{\underline{\ell}} = \hat{Y}_{\mathcal{J}, S}((t_{(\ell_1)}^{(S'_1)}, \dots, t_{(\ell_d)}^{(S'_d)})) $$

and because $t_{(\ell_1)}^{(S'_1)} = t_{(\ell_1)}^{(S_1)}$, we obtain from (36), $\hat{Y}_{\mathcal{J}, S}((t_{(\ell_1)}^{(S'_1)}, \dots, t_{(\ell_d)}^{(S'_d)})) = (\hat{\alpha}_{\mathcal{J}, S})_{\underline{\ell}}$ and thus

$$ \gamma_{\underline{\ell}} = (\hat{\alpha}_{\mathcal{J}, S})_{\underline{\ell}}. $$

Second, if $\ell_1 \in \{\nu+2, \dots, m_1+1\}$, we have

$$ \gamma_{\underline{\ell}} = \hat{Y}_{\mathcal{J}, S}((t_{(\ell_1)}^{(S'_1)}, \dots, t_{(\ell_d)}^{(S'_d)})) $$

and because $t_{(\ell_1)}^{(S'_1)} = t_{(\ell_1-1)}^{(S_1)}$, we obtain from (36),

$$ \hat{Y}_{\mathcal{J}, S}((t_{(\ell_1)}^{(S'_1)}, \dots, t_{(\ell_d)}^{(S'_d)})) = (\hat{\alpha}_{\mathcal{J}, S})_{(\ell_1-1, \ell_2, \dots, \ell_d)}. $$

Hence

$$ \gamma_{\underline{\ell}} = (\hat{\alpha}_{\mathcal{J}, S})_{(\ell_1-1, \ell_2, \dots, \ell_d)}. $$
Finally, if $\ell_1 = \nu + 1$, we have

$$ \gamma_{\underline{\ell}} = \hat{Y}_{\mathcal{J},S}((t_{(\nu+1)}^{(S'_1)}, t_{(\ell_2)}^{(S'_2)}, \dots, t_{(\ell_d)}^{(S'_d)})) \\ = (\hat{\alpha}_{\mathcal{J},S})(\nu, \ell_2, \dots, \ell_d) \phi_{t_{(\nu-1)}^{(S_1)}, t_{(\nu)}^{(S_1)}, t_{(\nu+1)}^{(S'_1)}}(t_{(\nu+1)}^{(S'_1)}) + (\hat{\alpha}_{\mathcal{J},S})(\nu+1, \ell_2, \dots, \ell_d) \phi_{t_{(\nu)}^{(S_1)}, t_{(\nu+1)}^{(S_1)}, t_{(\nu+2)}^{(S'_1)}}(t_{(\nu+1)}^{(S'_1)}), $$

from the definition of the multidimensional hat basis functions and the position of $(t_{(\nu+1)}^{(S'_1)}, t_{(\ell_2)}^{(S'_2)}, \dots, t_{(\ell_d)}^{(S'_d)})$ relatively to the knots in the subdivision $S$. Hence we have, since $t_{(\nu+1)}^{(S'_1)} = t$,

$$ \gamma_{\underline{\ell}} = (\hat{\alpha}_{\mathcal{J},S})(\nu, \ell_2, \dots, \ell_d) \frac{t_{(\nu+1)}^{(S_1)} - t}{t_{(\nu+1)}^{(S_1)} - t_{(\nu)}^{(S_1)}} + (\hat{\alpha}_{\mathcal{J},S})(\nu+1, \ell_2, \dots, \ell_d) \frac{t - t_{(\nu)}^{(S_1)}}{t_{(\nu+1)}^{(S_1)} - t_{(\nu)}^{(S_1)}}. $$

Hence, we have shown that

$$ \hat{Y}_{\mathcal{J},S} - \hat{Y}_{\mathcal{J},S'} = \sum_{\underline{\ell} \in \mathcal{L}_{S'}} \beta_{\underline{\ell}} \phi_{\underline{\ell}}^{(S')}, $$

with $\beta_{\underline{\ell}}$ as given in the proposition. To conclude the proof, by bilinearity of the square $L^2$ distance, it remains to prove that for $\underline{\ell} = (\ell_1, \dots, \ell_d) \in \mathcal{L}_{S'}$, $\underline{\ell}' = (\ell'_1, \dots, \ell'_d) \in \mathcal{L}_{S'}$,

$$ \int_{[0,1]^d} \phi_{\underline{\ell}}^{(S')}(x) \phi_{\underline{\ell}'}^{(S')}(x) dx = \Psi_{\underline{\ell}, \underline{\ell}'}, \quad (37) $$

with $\Psi_{\underline{\ell}, \underline{\ell}'}$ as given in the proposition. If there is $j \in \{1, \dots, d\}$ such that $|\ell_j - \ell'_j| > 1$, then the supports of $\phi_{\underline{\ell}}^{(S')}$ and $\phi_{\underline{\ell}'}^{(S')}$ are disjoint and thus (37) is zero. Consider now that for $j \in \{1, \dots, d\}$, $|\ell_j - \ell'_j| \le 1$. Then (37) is equal to

$$ \prod_{j=1, \dots, d} \int_0^1 \phi_{t_{(\ell_j-1)}^{(S'_j)}, t_{(\ell_j)}^{(S'_j)}, t_{(\ell_j+1)}^{(S'_j)}}(x) \phi_{t_{(\ell'_j-1)}^{(S'_j)}, t_{(\ell'_j)}^{(S'_j)}, t_{(\ell'_j+1)}^{(S'_j)}}(x) dx. $$

Let now $j \in \{1, \dots, d\}$ and write $t_{\ell}$ as a short-hand for $t_{(\ell)}^{(S'_j)}$, for $\ell \in \{0, \dots, m_{S'_j} + 1\}$. Write also $\phi_{\ell}$ as a short-hand for $\phi_{t_{(\ell-1)}^{(S'_j)}, t_{(\ell)}^{(S'_j)}, t_{(\ell+1)}^{(S'_j)}}$.

Consider $\ell, \ell' \in \{1, \dots, m_j + 1\}$. If $\ell = \ell' = 1$, we have

$$ \int_0^1 \phi_{\ell}(x) \phi_{\ell'}(x) = \int_{t_1}^{t_2} \left( \frac{t_2 - x}{t_2 - t_1} \right)^2 dx = \frac{t_2 - t_1}{3}. $$

If $\ell = \ell' = m_j + 1$, we have

$$ \int_0^1 \phi_{\ell}(x) \phi_{\ell'}(x) = \int_{t_{m_j}}^{t_{m_j+1}} \left( \frac{x - t_{m_j}}{t_{m_j+1} - t_{m_j}} \right)^2 dx = \frac{t_{m_j+1} - t_{m_j}}{3}. $$

If $\ell = \ell' \in \{2, \dots, m_j\}$, we have

$$ \int_0^1 \phi_{\ell}(x) \phi_{\ell'}(x) = \int_{t_{\ell-1}}^{t_{\ell}} \left( \frac{x - t_{\ell-1}}{t_{\ell} - t_{\ell-1}} \right)^2 dx + \int_{t_{\ell}}^{t_{\ell+1}} \left( \frac{t_{\ell+1} - x}{t_{\ell+1} - t_{\ell}} \right)^2 dx = \frac{t_{\ell+1} - t_{\ell-1}}{3}. $$

If $\ell' = \ell + 1 \in \{2, \dots, m_j + 1\}$, we have

$$ \int_0^1 \phi_{\ell}(x) \phi_{\ell'}(x) = \int_{t_{\ell}}^{t_{\ell+1}} \left( \frac{t_{\ell+1} - x}{t_{\ell+1} - t_{\ell}} \right) \left( \frac{x - t_{\ell}}{t_{\ell+1} - t_{\ell}} \right) dx = \frac{t_{\ell+1} - t_{\ell}}{6}. $$

Finally, if $\ell' = \ell - 1 \in \{1, \dots, m_j\}$, we have

$$ \int_0^1 \phi_{\ell}(x) \phi_{\ell'}(x) = \int_{t_{\ell-1}}^{t_{\ell}} \left( \frac{t_{\ell} - x}{t_{\ell} - t_{\ell-1}} \right) \left( \frac{x - t_{\ell-1}}{t_{\ell} - t_{\ell-1}} \right) dx = \frac{t_{\ell} - t_{\ell-1}}{6}. $$

Hence, we have shown that for $j \in \{1, \dots, d\}$ and $\ell_j, \ell'_j \in \{1, \dots, m_j + 1\}$ with $|\ell_j - \ell'_j| \le 1$,

$$ \int_0^1 \phi_{t_{(\ell_j-1)}^{(S'_j)}, t_{(\ell_j)}^{(S'_j)}, t_{(\ell_j+1)}^{(S'_j)}}(x) \phi_{t_{(\ell'_j-1)}^{(S'_j)}, t_{(\ell'_j)}^{(S'_j)}, t_{(\ell'_j+1)}^{(S'_j)}}(x) dx = \psi_{\ell_j, \ell'_j}^{(j)}. $$

with the notation of the proposition. This concludes the proof. $\square$
**Proof of Proposition 4.** The current mode function

$$ \hat{Y}_{\mathcal{J}, S} = \sum_{\tilde{\ell} \in \mathcal{L}_S} (\hat{\alpha}_{\mathcal{J}, S})_{\tilde{\ell}} \phi_{\tilde{\ell}}^{(S)} $$

can be expressed of the form

$$ \hat{Y}_{\mathcal{J}, S} = \sum_{\tilde{\ell} \in \mathcal{L}_S} \sum_{\ell_{d+1}=1}^{2} (\hat{\alpha}_{\mathcal{J}, S})_{\tilde{\ell}} \phi_{\tilde{\ell}}^{(S)} \phi_{t_{(\ell_{d+1}-1)}, t_{(\ell_{d+1})}, t_{(\ell_{d+1})}}^{(S^0)} = \sum_{\ell \in \mathcal{L}_{S'}} \gamma_{\ell} \phi_{\ell}^{(S+(d+1))} $$

with $\gamma_{\ell} = (\hat{\alpha}_{\mathcal{J}, S})_{\tilde{\ell}}$ when $\ell = (\tilde{\ell}, \ell_{d+1})$. Hence we obtain

$$ \hat{Y}_{\mathcal{J}, S} - \hat{Y}_{\mathcal{J} \cup \{d+1\}, S+(d+1)} = \sum_{\ell \in \mathcal{L}_{S+(d+1)}} \beta_{\ell} \phi_{\ell}^{(S+(d+1))}, $$

with $\beta_{\ell}$ as defined in the proposition. The rest of the proof is then identical to the proof of Proposition 3. $\square$

**Proof of Proposition 1.** Let us first prove that

$$ L_1 \dots L_d f(t) = \sum_{\epsilon_1, \dots, \epsilon_d \in \{-, +\}} \left( \prod_{j=1}^{d} \omega_{\epsilon_j}(t_j) \right) f(t_1^{\epsilon_1}, \dots, t_d^{\epsilon_d}). $$

In dimension 1, notice that the affine extension can be rewritten as

$$ L_B f(t) = \omega_-(t) f(t^-) + \omega_+(t) f(t^+) = \sum_{\epsilon \in \{-, +\}} \omega_\epsilon(t) f(t^\epsilon), \quad (38) $$

with the convention chosen for $\omega_-(t)$, $\omega_+(t)$ for $t \in B$. Indeed, in that case $t^- = t^+ = t$, and thus $L_B f(t) = \frac{1}{2}(f(t) + f(t)) = f(t)$. Then, let us show by induction on $i = d, d-1, \dots, 1$ the property

$$ (P_i): \quad \forall t \in F^{1:i-1} \times [0, 1]^{d-i+1}, \\ L_i \dots L_d f(t) = \sum_{\epsilon_i, \dots, \epsilon_d \in \{-, +\}} \left( \prod_{j=i}^{d} \omega_{\epsilon_j}(t_j) \right) f(t_1, \dots, t_{i-1}, t_i^{\epsilon_i}, \dots, t_d^{\epsilon_d}). \quad (39) $$

First $P_d$ is true, by the expression (38) for $L_{F_d}$. Now, assume that $P_{i+1}$ is true (for $i \in \{1, \dots, d-1\}$). Then, we have, for $t \in F^{1:i-1} \times [0, 1]^{d-i+1}$:

$$ \begin{aligned} L_i \dots L_d f(t) &= L_{F_i} (L_{i+1} \dots L_d f(\cdot, t_{\sim i}))(t_i) \\ &= \sum_{\epsilon_i \in \{-, +\}} \omega_{\epsilon_i}(t_i) (L_{i+1} \dots L_d f)(t_i^{\epsilon_i}, t_{\sim i}) \\ &= \sum_{\epsilon_i \in \{-, +\}} \omega_{\epsilon_i}(t_i) \left( \sum_{\epsilon_{i+1}, \dots, \epsilon_d \in \{-, +\}} \left( \prod_{j=i+1}^{d} \omega_{\epsilon_j}(t_j) \right) f(t_1, \dots, t_{i-1}, t_i^{\epsilon_i}, t_{i+1}^{\epsilon_{i+1}}, \dots, t_d^{\epsilon_d}) \right) \end{aligned} $$

which gives $P_i$. Finally $P_1$ gives (18). Furthermore, notice that:

$$ \sum_{\epsilon_1, \dots, \epsilon_d \in \{-, +\}} \left( \prod_{j=1}^{d} \omega_{\epsilon_j}(t_j) \right) = \prod_{j=1}^{d} (\omega_-(t_j) + \omega_+(t_j)) = 1. $$

Let us now justify the existence and unicity of $P_{F \to [0,1]^d}$.

Firstly, if $g$ exists, then necessarily $g = L_1 \dots L_d f$. Indeed, by assumption, when $t_{\sim 1} = (t_2, \dots, t_d)$ is fixed in $[0, 1]^{d-1}$, the univariate function

$$ g(\cdot, t_{\sim 1}) : u_1 \mapsto g(u_1, t_{\sim 1}) $$

is continuous and affine on all intervals of $[0, 1] \setminus F_1$. By property of the affine extension, it is equal to the affine extension of its restriction to $F_1$:

$$ g(\cdot, t_{\sim 1}) = L_{F_1} (g(\cdot, t_{\sim 1})|_{F_1}) = L_{F_1} (g|_{F_1 \times [0,1]^{d-1}}(\cdot, t_{\sim 1})). $$
This shows that

$$g = L_1 g_{|F_1 \times [0,1]^{d-1}}.$$

By the same reasoning,

$$g_{|F_1 \times [0,1]^{d-1}} = L_2 g_{|F_1 \times F_2 \times [0,1]^{d-2}},$$

and by an immediate induction,

$$g = L_1 \dots L_d g_F = L_1 \dots L_d f.$$

Now, let us check that the function $g = L_1 \dots L_d f$ verifies the conditions of the proposition.

* It is equal to $f$ on $F$, since each $L_i$ leaves the values of its input function unchanged on $F$.

* Let check that all 1-dimensional cuts of $g$ are affine on all intervals of the complements of the $F_i$'s. Indeed, consider the explicit formula (18). Let fix $i \in \{1, \dots, d\}$ and consider an interval $I_i$ included in $[0, 1] \setminus F_i$. Without loss of generality, we assume that $I_i = [a_i, b_i]$ is closed, with $a_i < b_i$. Then, for all $t_i \in I_i$, we have $a_i^- = t_i^- < t_i^+ = b_i^+$. Thus, $t_i^-, t_i^+$ do not depend on $t_i$. Consequently, $\omega_+(t_i) = \frac{t_i - a_i^-}{b_i^+ - a_i^-}$ and $\omega_-(t_i) = 1 - \omega_+(t_i)$ depend linearly on $t_i$. Finally, by (18), for all $t_{i-} \in [0, 1]^{d-1}$, the cut function $t_i \mapsto L_1 \dots L_d f(t)$ is affine on $I_i$.

* To prove continuity, by composition, it is sufficient to prove that for $i \in \{1, \dots, d\}$, when $h$ belongs to $C(F^{1:i} \times [0, 1]^{d-i}, \mathbb{R})$ then $L_i h$ belongs to $C(F^{1:(i-1)} \times [0, 1]^{d-i+1}, \mathbb{R})$. Let us thus consider $i$, and $h$ as just described. Let $t = (t_1, \dots, t_d) \in F^{1:(i-1)} \times [0, 1]^{d-i+1}$ and consider a sequence $t_n = (t_{n,1}, \dots, t_{n,d}) \in F^{1:(i-1)} \times [0, 1]^{d-i+1}$ converging to $t$. We will show that $L_i h(t_n)$ converges to $L_i h(t)$.

Up to extracting subsequences, it suffices to consider the cases (1) $t_{n,i} \in F_i$, (2) $t_{n,i} \notin F_i$, $t_{n,i} \le t_i$, $t_{n,i}$ increasing and (3) $t_{n,i} \notin F_i$, $t_{n,i} \ge t_i$, $t_{n,i}$ decreasing.

In case (1), since $F_i$ is closed we have $t_i \in F_i$. Thus, as $h$ is continuous, $L_i h(t_n) = h(t_n)$ converges to $h(t) = L_i h(t)$.

Consider now the case (2). Observe that we have

$$L_i h(t_n) = [1 - \omega_+(t_{n,i})]h(t_{n,i}^-, t_{n,i-}) + \omega_+(t_{n,i})h(t_{n,i}^+, t_{n,i-}). \quad (40)$$

Furthermore, as $t_{n,i}^-$ is increasing and bounded above by $t_i$, it converges to a limit $t_{\infty,i}^- \le t_i$.

Consider the case (2a) where $t_i \in F_i$. Then, by definition of $t_{n,i}^+$, we have $t_{n,i} \le t_{n,i}^+ \le t_i$.

Hence $t_{n,i}^+$ converges to $t_i$. If, first, $t_{\infty,i}^- = t_i$, then from (40), $L_i h(t_n)$ is a convex combination of two values of $h$ at two inputs that converge to $t$ so $L_i h(t_n)$ converges to $h(t) = L_i h(t)$:

$$\begin{aligned} & |L_i h(t_n) - h(t)| \\ &= |[1 - \omega_+(t_{n,i})](h(t_{n,i}^-, t_{n,i-}) - h(t)) + \omega_+(t_{n,i})(h(t_{n,i}^+, t_{n,i-}) - h(t))| \\ &\le \max(|h(t_{n,i}^-, t_{n,i-}) - h(t)|, |h(t_{n,i}^+, t_{n,i-}) - h(t)|) \underset{n \to +\infty}{\longrightarrow} 0. \end{aligned}$$

If, second, $t_{\infty,i}^- < t_i$, then $\omega_+(t_{n,i}) \to 1$ and by (40), $L_i h(t_n)$ converges to $h(t) = L_i h(t)$.

Consider now the case (2b) where $t_i \notin F_i$. Then, since $[0, 1] \setminus F_i$ is open in $[0, 1]$, for $n$ large enough we have $t_{n,i}^- = t_i^-$ and $t_{n,i}^+ = t_i^+$. Thus, $\omega_+(t_{n,i}) \to \omega_+(t_i)$, and by (40), $L_i h(t_n) \to L_i h(t)$. This concludes the case (2). The case (3) is treated similarly.

□

### Proof of proposition 2.

1. First, the linearity of $P_{F \to [0,1]^d}$ comes from the linearity of $g \mapsto L_B g$, and by composition of affine maps.

Now, let us prove that $P_{F \to [0,1]^d}$ is 1-Lipschitz. By linearity, it is sufficient to show that for $f \in C(F, \mathbb{R})$,

$$\sup_{t \in [0,1]^d} |P_{F \to [0,1]^d}(f)(t)| \le \sup_{t \in F} |f(t)|. \quad (41)$$

Notice that, with $g$ a univariate and continuous function on a closed subset $B$ of $[0, 1]$ containing the boundaries 0 and 1, by construction, the values of the affine extension of $L_B(g)$ lie in the range of $g$ values. This implies that, for $f \in C(F, \mathbb{R})$,

$$\sup_{t \in F_1 \times \cdots \times F_{d-1} \times [0,1]} |L_d f(t)| \le \sup_{t_{\sim d} \in F_1 \times \cdots \times F_{d-1} \atop u_d \in F_d} |f(u_d, t_{\sim d})| \le \sup_{t \in F} |f(t)|.$$
Similarly, we have

$$
\sup_{t \in F_1 \times \cdots \times F_{d-2} \times [0,1]^2} |L_{d-1} L_d f(t)| \le \sup_{\substack{t \sim (d-1) \in F_1 \times \cdots \times F_{d-2} \times [0,1] \\ u_{d-1} \in F_{d-1}}} |L_d f(u_{d-1}, t \sim (d-1))| \\
\le \sup_{t \in F} |f(t)|.
$$

Hence, by iteration we show (41).

2. Any $f \in E_S$ can be written as $Y_{S,\alpha}$ with $\alpha \in A_S$. Let us first consider the marginal extension of $Y_{S,\alpha|F}$ with respect to coordinate $d$: for $t = (t_1, \dots, t_d) \in F_1 \times \dots \times F_{d-1} \times [0, 1]$, we have

$$
L_d Y_{S,\alpha|F}(t) = L_{F_d}[u_d \mapsto Y_{S,\alpha|F}(t_1, \dots, t_{d-1}, u_d)](t_d).
$$

Observe that, when $t_1, \dots, t_{d-1}$ are fixed in $F_1 \times \dots \times F_{d-1}$, the univariate function

$$
u_d \mapsto Y_{S,\alpha|F_1 \times \dots \times F_{d-1} \times [0,1]}(t_1, \dots, t_{d-1}, u_d)
$$

is piecewise linear, and all its knots are contained in $F_d$, by definition of $F$. Hence, when restricting it to $F_d$ and then reextending it to $[0, 1]$ with $L_{F_d}$, one obtains exactly the same function. In other words: for $t = (t_1, \dots, t_d) \in F_1 \times \dots \times F_{d-1} \times [0, 1]$, we have

$$
L_d Y_{S,\alpha|F}(t) = Y_{S,\alpha|F_1 \times \dots \times F_{d-1} \times [0,1]}(t).
$$

By the same reasoning, we have, for $t = (t_1, \dots, t_d) \in F_1 \times \dots \times F_{d-2} \times [0, 1]^2$,

$$
L_{d-1} L_d Y_{S,\alpha|F}(t) = Y_{S,\alpha|F_1 \times \dots \times F_{d-2} \times [0,1]^2}(t).
$$

Hence by an immediate induction: $P_{F \to [0,1]^d} Y_{S,\alpha|F}(t) = Y_{S,\alpha}(t)$ for $t \in [0,1]^d$.

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

**Proof of Corollary 1.** The functions $f$ and $g$ are $d$-affine on $\Delta$, thus they are polynomial functions and can be defined on $[0, 1]^d$. Consider $F = \prod_{j=1}^d ([0, 1] \setminus (x_j^-, x_j^+))$. One can simply show that $f$ satisfies the conditions for $P_{F \to [0,1]^d}(f)$ in Proposition 1, then by unicity $f = P_{F \to [0,1]^d}(f)$. Similarly $g = P_{F \to [0,1]^d}(g)$. Finally, from (18), $P_{F \to [0,1]^d}(g) = P_{F \to [0,1]^d}(f)$ since $f$ and $g$ coincide on the $2^d$ vertices of $\Delta$. $\square$

**Proof of Lemma 2.** Observe that $\dot{C} = \text{int}_{|||\infty}(C)$. Let us consider $g \in \mathcal{H} \cap \text{int}_{|||\infty}(C)$. Then $g \in \mathcal{H} \cap C$. Since $k$ is continuous and defined on the compact set $[0,1]^d \times [0,1]^d$, there exists a constant $C_{\text{sup}}$ such that for all $h \in C([0,1]^d, \mathbb{R})$, $||h||_{\infty} \le C_{\text{sup}}||h||_{\mathcal{H}}$ (see for instance Lemma 2 in [5]). Let $\epsilon > 0$ such that $||g-h||_{\infty} \le \epsilon$ implies $h \in C$. Then, for all $h \in \mathcal{H}$ such that $||g-h||_{\mathcal{H}} \le \epsilon/C_{\text{sup}}$, we have $h \in C$. This means that $g$ is in $\text{int}_{|||\infty}(\mathcal{H} \cap C)$. $\square$

**Lemma 4.** Let $B$ be a closed subset of $[0, 1]$ containing 0 and 1. Let $f_B \in C(B, \mathbb{R})$ be non-decreasing. Then $L_B f_B$ is non-decreasing from $[0, 1] \to \mathbb{R}$.

**Proof of Lemma 4.** Write $L = L_B f_B$. Let $0 \le v < w \le 1$. Since $L$ is affine on $[v^-, v^+]$ and $L(v^-) = f_B(v^-) \le f_B(v^+) = L(v^+)$, we have $L(v^-) \le L(v) \le L(v^+)$. Similarly $L(w^-) \le L(w) \le L(w^+)$. Hence if (1) $v^+ \le w^-$, then $L(v) \le L(w)$. If (2) $v^+ > w^-$, then one can show that $[v, w] \cap B = \emptyset$ and thus $v^- = w^-$ and $v^+ = w^+$. Then $L$ is affine on $[v^-, v^+]$ with $L(v^-) \le L(v^+)$ and $v, w \in [v^-, v^+]$. Hence $L(v) \le L(w)$. Hence $L$ is increasing from $[0, 1] \to \mathbb{R}$. $\square$

In the next lemma, for $B \subseteq \mathbb{R}$, we assume that $f_B : B \to \mathbb{R}$ is convex, i.e. for all $x, y \in B$, $\lambda \in [0, 1]$ with $\lambda x + (1 - \lambda)y \in B$, we have

$$
f_B(\lambda x + (1 - \lambda)y) \le \lambda f_B(x) + (1 - \lambda)f_B(y). \quad (42)
$$

**Lemma 5.** Let $B$ be a closed subset of $[0, 1]$ containing 0 and 1. Let $f_B \in C(B, \mathbb{R})$ be convex. Then $L_B f_B$ is convex from $[0, 1] \to \mathbb{R}$.

**Proof of Lemma 5.** For a function $g$ defined on a domain $D \subseteq \mathbb{R}$, we recall the definition of the epigraph of $f$, denoted $\text{epi}(f)$:

$$
\text{epi}(f) := \{(t, y) \in D \times \mathbb{R} \text{ s.t. } y \ge f(t)\}.
$$
We also recall that the convex hull of a set $A \subseteq \mathbb{R}^2$, denoted hull(A), is the smallest convex set containing A. Equivalently, it is the set of convex combinations of points in A (see e.g. [22]). Then, we are going to prove that, if $f_B$ is convex, then the convex hull of the epigraph of $f_B$ is the epigraph of its affine extension $L_B f_B$:

$$ \text{hull}(\text{epi}(f_B)) = \text{epi}(L_B f_B). $$

The result will then be deduced since it shows that $\text{epi}(L_B f_B)$ is convex, which is equivalent to the convexity of $L_B f_B$ (as $L_B f_B$ is defined on the interval $[0, 1]$). Let us now come back to the proof.

* First let us show that (and even if $f_B$ is not convex) $\text{epi}(L_B f_B) \subseteq \text{hull}(\text{epi}(f_B))$. Let $(t, y) \in \text{epi}(L_B f_B)$.

- If $t \in B$, then $y \ge L_B f_B(t) = f_B(t)$. Thus $(t, y) \in \text{epi}(f_B) \subseteq \text{hull}(\text{epi}(f_B))$.

- If $t \notin B$, then $t^- < t < t^+$. Consider the straight line joining $(t^-, f_B(t^-))$ and $(t, y)$. Let $(t^+, y^+)$ the point on that straight line with abscissa $t^+$. Notice that $L_B f_B$ is a straight line on $[t^-, t^+]$, joining $(t^-, f_B(t^-))$ and $(t^+, f_B(t^+))$. Then, by Thalès theorem, the sign of $y^+ - f_B(t^+)$ is the same as $y - L_B f_B(t)$, which is positive because $(t, y) \in \text{epi}(L_B f_B)$. Thus $(t^+, y^+) \in \text{epi}(f_B)$. Finally, we have shown that $(t, y)$ belongs to a line segment whose endpoints $(t^-, f_B(t^-))$, $(t^+, y^+)$ belong to $\text{epi}(f_B)$, which proves that $(t, y) \in \text{hull}(\text{epi}(f_B))$.

* Conversely, let us prove that $\text{hull}(\text{epi}(f_B)) \subseteq \text{epi}(L_B f_B)$. Here the convexity of $f_B$ is required.

Let $(t, y) \in \text{hull}(\text{epi}(f_B))$. Thus $(t, y)$ belongs to a polygon $P$, which is either a singleton, a segment or a triangle and whose vertices $(t_1, y_1), \dots, (t_m, y_m)$ are in $\text{epi}(f_B)$, with $m \in \{1, 2, 3\}$, from Carathéodory's theorem. The intersection of $P$ with the band $[t^-, t^+] \times \mathbb{R}$ is a convex polygon containing $(t, y)$. Thus, from Krein-Milman theorem, $(t, y)$ is a convex combination of its extremal points $(u_1, z_1), \dots, (u_r, z_r)$. By definition of $t^-, t^+$, the points $t_1, \dots, t_m \notin (t^-, t^+)$. Then, one can see that these extremal points have abscissas in $\{t^-, t^+\}$ and either belong to $\{(t_1, y_1), \dots, (t_m, y_m)\}$ or are in the segments with endpoints in $\{(t_1, y_1), \dots, (t_m, y_m)\}$.

Let us now prove that these extremal points are in $\text{epi}(f_B)$. Consider for instance $(u_1, z_1)$.

- If $(u_1, z_1)$ is in $\{(t_1, y_1), \dots, (t_m, y_m)\}$, it belongs to $\text{epi}(f_B)$.

- Otherwise, assume without loss of generality that $(u_1, z_1)$ is the intersection of the segment joining $(t_1, y_1)$, $(t_2, y_2)$ and $\{t^-\} \times \mathbb{R}$, with $t_1 \le t^- \le t_2$. Then, there exists $\lambda \in [0, 1]$ such that $(u_1, z_1) = \lambda(t_1, y_1) + (1 - \lambda)(t_2, y_2)$. We have, by using that $(t_1, y_2), (t_2, y_2) \in \text{epi}(f_B)$ and by convexity of $f_B$,

$$ \begin{aligned} z_1 &= \lambda y_1 + (1 - \lambda) y_2 \\ &\ge \lambda f_B(t_1) + (1 - \lambda) f_B(t_2) \\ &\ge f_B(\lambda t_1 + (1 - \lambda) t_2) \\ &= f_B(u_1). \end{aligned} $$

Hence $(u_1, z_1)$ belongs to $\text{epi}(f_B)$.

Now, $L_B f_B$ is affine on $[t^-, t^+]$. Recall that $u_1, \dots, u_r$ are elements of $\{t^-, t^+\}$, and thus are in $B$. Hence $L_B f_B$ coincides with $f_B$ at $u_1, \dots, u_m$. Writing that $(t, y)$ is a convex combination of $(u_1, z_1), \dots, (u_r, z_r)$, and using that $(u_1, z_1), \dots, (u_r, z_r) \in \text{epi}(f_B)$, we obtain:

$$ y = \sum_{i=1}^{r} \lambda_i z_i \ge \sum_{i=1}^{r} \lambda_i f_B(u_i) = \sum_{i=1}^{r} \lambda_i L_B f_B(u_i) = L_B f_B \left( \sum_{i=1}^{r} \lambda_i u_i \right) = L_B f_B(t). $$

This proves that $(t, y) \in \text{epi}(L_B f_B)$.

□

**Proof of Lemma 3.** Consider that $\mathcal{C}$ is given by (14). Let $f \in \mathcal{C}$, with $f|_F$ the restriction of $f$ to $F$. From (18), for $x = (x_1, \dots, x_d) \in [0, 1]^d$, the value of $P_{F \to [0,1]^d}(f|_F)(x)$ is a convex combination of the values $(f(x_1^{\epsilon_1}, \dots, x_d^{\epsilon_d}))_{\epsilon_1, \dots, \epsilon_d \in \{-, +\}}$. These $2^d$ values are in $[a, b]$ and thus $P_{F \to [0,1]^d}(f|_F)(x)$ is also in $[a, b]$. Thus $P_{F \to [0,1]^d}(f|_F) \in \mathcal{C}$.
Consider that $\mathcal{C}$ is given by (15). Let $i \in \{1, \dots, d\}$. Let $x \in [0, 1]^d$ and $0 \le u_i < v_i \le 1$. From (18), we have, with the notation $x_{\sim i}^{\epsilon} = (x_1^{\epsilon_1}, \dots, x_{i-1}^{\epsilon_{i-1}}, x_{i+1}^{\epsilon_{i+1}}, \dots, x_d^{\epsilon_d})$,

$$
\begin{align}
& P_{F \to [0,1]^d}(f_F)(v_i, x_{\sim i}) - P_{F \to [0,1]^d}(f_F)(u_i, x_{\sim i}) \nonumber \\
&= \sum_{\epsilon_1, \dots, \epsilon_d \in \{-,+\}} \left( w_{\epsilon_i}(v_i) \prod_{\substack{j=1 \\ j \neq i}}^d \omega_{\epsilon_j}(x_j) \right) f(v_i^{\epsilon_i}, x_{\sim i}^{\epsilon}) \nonumber \\
&\quad - \sum_{\epsilon_1, \dots, \epsilon_d \in \{-,+\}} \left( w_{\epsilon_i}(u_i) \prod_{\substack{j=1 \\ j \neq i}}^d \omega_{\epsilon_j}(x_j) \right) f(u_i^{\epsilon_i}, x_{\sim i}^{\epsilon}) \nonumber \\
&= \sum_{\epsilon_{\sim i} \in \{-,+\}^{d-1}} \left( \prod_{\substack{j=1 \\ j \neq i}}^d \omega_{\epsilon_j}(x_j) \left( \sum_{\epsilon_i \in \{-,+\}} w_{\epsilon_i}(v_i) f(v_i^{\epsilon_i}, x_{\sim i}^{\epsilon}) - \sum_{\epsilon_i \in \{-,+\}} w_{\epsilon_i}(u_i) f(u_i^{\epsilon_i}, x_{\sim i}^{\epsilon}) \right) \right) \nonumber \\
&= \sum_{\epsilon_{\sim i} \in \{-,+\}^{d-1}} \left( \prod_{\substack{j=1 \\ j \neq i}}^d \omega_{\epsilon_j}(x_j) \left( L_{F_i} f_{|F}(\cdot, x_{\sim i}^{\epsilon})(v_i) - L_{F_i} f_{|F}(\cdot, x_{\sim i}^{\epsilon})(u_i) \right) \right). \tag{43}
\end{align}
$$

In the above display, the function $f_{|F}(\cdot, x_{\sim i}^{\epsilon})$ is continuous non-decreasing from $F_i$ to $\mathbb{R}$. Hence, from Lemma 4, the function

$t \in [0, 1] \mapsto L_{F_i} f_{|F}(\cdot, x_{\sim i}^{\epsilon})(t)$

is non-decreasing. Hence, the difference in (43) is non-negative and the weighted sum is non-negative since the weights are non-negative. This concludes the proof.

Consider that $\mathcal{C}$ is given by (16). Then the proof is identical to the case where $\mathcal{C}$ is given by (15). Instead of using Lemma 4 we use Lemma 5. For $i \in \{1, \dots, d\}$, $x \in [0, 1]^d$ and $0 \le u_i < v_i < w_i \le 1$, instead of considering the first order finite difference

$$P_{F \to [0,1]^d}(f_F)(v_i, x_{\sim i}) - P_{F \to [0,1]^d}(f_F)(u_i, x_{\sim i}),$$

we consider the second order finite difference

$$
\frac{P_{F \to [0,1]^d}(f_F)(w_i, x_{\sim i}) - P_{F \to [0,1]^d}(f_F)(v_i, x_{\sim i})}{w_i - v_i} - \frac{P_{F \to [0,1]^d}(f_F)(v_i, x_{\sim i}) - P_{F \to [0,1]^d}(f_F)(u_i, x_{\sim i})}{v_i - u_i}.
$$

□

**Proof of Theorem 2.** To prove the uniform convergence, we cannot directly apply Theorem 1 because the set of d-dimensional knots is not assumed to be dense in $[0, 1]^d$. The idea is to match this situation, by considering functions defined on their closure $F$ and by using the multiaffine extension $P_{F \to [0,1]^d}$.

We will thus find functions and sets corresponding to $\hat{Y}_{S^{(m)}}$, $k$, $\mathcal{H}$, $\mathcal{C}$ and $\mathcal{I}_{X,y^{(n)}}$ in Theorem 1. We will then show that conditions corresponding to those of Theorem 1 hold. Under these new conditions, the proof of Theorem 1 can be repeated, enabling us to obtain the conclusion of Theorem 2.

To $\hat{Y}_{S^{(m)}}$, we thus associate $\hat{Y}_{S^{(m)}|F}$, its restriction to $F$. The input space of $\hat{Y}_{S^{(m)}}$ is $[0, 1]^d$ (written $X$ in [6]) while the input space of $\hat{Y}_{S^{(m)}|F}$ is $F$. To $k$ we associate $k_F$. To $\mathcal{I}_{X,y^{(n)}}$ (the set of equality constraints, written $I$ in [6]) we associate $\mathcal{I}_{F,X,y^{(n)}}$. To $\mathcal{C}$ (the set of inequality constraints) we associate $\mathcal{C}_F$. To $\phi_{\underline{\ell}}^{(S^{(m)})}$, $\underline{\ell} \in \mathcal{L}_{S^{(m)}}$ (the basis functions) we associate their restriction to $F$, $\phi_{\underline{\ell}|F}^{(S^{(m)})}$. Then notice that $\hat{Y}_{S^{(m)}|F}$ indeed corresponds to the mode function $\hat{Y}_{S^{(m)}}$ of Theorem 1, in the sense that we have, for $t \in F$,

$$\hat{Y}_{S^{(m)}|F} = \sum_{\underline{\ell} \in \mathcal{L}_{S^{(m)}}} (\hat{\alpha}_{S^{(m)}})_{\underline{\ell}} \phi_{\underline{\ell}|F}^{(S^{(m)})}(t)$$
and

$$
\begin{align}
\hat{\alpha}_{S^{(m)}} &\in \underset{\alpha \in A_{S^{(m)}}}{\operatorname{argmin}} \quad \alpha^{\top} k(S^{(m)}, S^{(m)})^{-1}\alpha \nonumber \\
&\qquad Y_{S^{(m)}, \alpha} \in \mathcal{I}_{X,y(n)} \cap \mathcal{C} \nonumber \\
&= \underset{\alpha \in A_{S^{(m)}}}{\operatorname{argmin}} \quad \alpha^{\top} k_F(S^{(m)}, S^{(m)})^{-1}\alpha. \tag{44}
\end{align}
$$

Above, we have used that $P_{F \to [0,1]^d} Y_{S^{(m)}, \alpha|F} = Y_{S^{(m)}, \alpha}$, from Proposition 2. We finally remark that, by construction, $F$ contains all the knots $t_{\ell}^{(S^{(m)})}$, $\ell \in \mathcal{L}_{S^{(m)}}$ obtained from $S^{(m)}$.

Then, let us prove that the following assumptions, corresponding to those required for Theorem 1, are satisfied. Notice that $\mathcal{C}_F$ is convex, since the map $P_{F \to [0,1]^d}$ is affine and $\mathcal{C}$ is convex. Then, we will show that the following extensions to (H1) and (H2) hold:

$$
(H1, F) \int_{||.||_{\mathcal{H}_F}} (\mathcal{H}_F \cap \mathcal{C}_F) \cap \mathcal{I}_{F,X,y(n)} \neq \emptyset.
$$

$$
(H2, F) \quad \pi_{S^{(m)}|F}(\mathcal{C}_F) \subseteq \mathcal{C}_F
$$

Let us first prove (H2, F). Let $f_F \in \mathcal{C}_F$. Thus $f = P_{F \to [0,1]^d}(f_F) \in \mathcal{C}$. Therefore $\pi_{S^{(m)}}(f) \in \mathcal{C}$ from Condition 5. Now $\pi_{S^{(m)}|F}(f_F)$ is the restriction to $F$ of the function $\pi_{S^{(m)}}(f)$, since $f = f_F$ on $F$ (by definition of $P_{F \to [0,1]^d}$). Finally, by Proposition 2,

$$
P_{F \to [0,1]^d} (\pi_{S^{(m)}|F}(f_F)) = P_{F \to [0,1]^d} ((\pi_{S^{(m)}}(f))_{|F}) = \pi_{S^{(m)}}(f) \in \mathcal{C}.
$$

Hence $\pi_{S^{(m)}|F}(f_F) \in \mathcal{C}_F$ which shows (H2, F).

Let us now prove (H1, F). For $i = 1, \dots, n$ let us write $x^{(i)} = (x_1^{(i)}, \dots, x_d^{(i)})$. For $j = 1, \dots, d$, write $x_{j,-}^{(i)} = \max\{u \in F_j; u \le x_j^{(i)}\}$ and $x_{j,+} = \min\{u \in F_j; u \ge x_j^{(i)}\}$. The set

$$
\{(x_{1,\epsilon_1}^{(i)}, \dots, x_{d,\epsilon_d}^{(i)})\}_{\epsilon_1, \dots, \epsilon_d \in \{-+,\}, i=1,\dots,n}
$$

can be written as $\{w_1, \dots, w_p\}$ with $w_1, \dots, w_p$ two by two distinct.

From Condition 2, for each $m \ge 0$ we have

$$
\emptyset \neq \{\alpha \in A_{S^{(m_0)}}, Y_{S^{(m_0)}, \alpha} \in \mathcal{I}_{X,y(n)} \cap \mathcal{C}^{\bullet}\} \subseteq \{\alpha \in A_{S^{(m)}}, Y_{S^{(m)}, \alpha} \in \mathcal{I}_{X,y(n)} \cap \mathcal{C}^{\bullet}\}, \quad (45)
$$

because the sequence of function spaces ${Y_{S^{(m)},\alpha}; \alpha \in A_{S^{(m)}}}_{m \ge m_0}$ is nested. Hence, we can take $\alpha$ in the set in (45). Write, for $i = 1, \dots, p$, $z_i = Y_{S^{(m)},\alpha}(w_i)$.

Then from Condition 4, the set $\text{int}_{||.||_{\mathcal{H}}}(\mathcal{H} \cap \mathcal{C}) \cap \mathcal{I}_{W,z(n)}$ is non-empty. For $f$ in this set, let us show that $f|_F$ belongs to $\text{int}_{||.||_{\mathcal{H}_F}}(\mathcal{H}_F \cap \mathcal{C}_F) \cap \mathcal{I}_{F,X,y(n)}$. Thus, this set will be non-empty, what is to prove.

Firstly, let us check that $f|_F \in \mathcal{I}_{F,X,y(n)}$. Indeed, on each of the $n$ hypercubes

$$
\prod_{j=1,...,d} [x_{j,-}^{(i)}, x_{j,+}^{(i)}],
$$

$i = 1, \dots, n$, $P_{F \to [0,1]^d}(f_F)$ coincide with $Y_{S^{(m)},\alpha}$. Indeed, consider one of these hypercubes. The 2$^d$ vertices of this hypercube belong to $F$, so on these 2$^d$ points, $P_{F \to [0,1]^d}(f_F)$ coincide with $f|_F$ which coincide with $f$ which coincide with $Y_{S^{(m)},\alpha}$. Furthermore, the two functions $P_{F \to [0,1]^d}(f_F)$ and $Y_{S^{(m)},\alpha}$ are $d$-affine on this hypercube, and we have shown that they take the same values on the 2$^d$ vertices. Thus they are equal on the hypercube from Corollary 1. Hence in particular $P_{F \to [0,1]^d}(f_F)(x^{(i)}) = Y_{S^{(m)},\alpha}(x^{(i)}) = y_i$ and thus $P_{F \to [0,1]^d}(f_F) \in \mathcal{I}_{X,y(n)}$ and so $f|_F \in \mathcal{I}_{F,X,y(n)}$. Note that the above argumentation still goes through in the case where there exist $i, j$'s such that $x_{j,-}^{(i)} = x_{j,+}^{(i)}$.

Secondly, notice that $f \in \mathcal{H} \cap \mathcal{C}$. By Theorem 6 in [10], $\mathcal{H}_F$ is formed by restrictions to $F$ of functions in $\mathcal{H}$, and for all $f_F \in \mathcal{H}_F$, we have $\|f_F\|_{\mathcal{H}_F} = \inf_{h \in \mathcal{H}, h|_F = f_F} \|h\|$. Thus $f|_F \in \mathcal{H}_F$. Moreover, $f|_F \in \mathcal{C}_F$ by Condition 6. Furthermore, let $\epsilon > 0$ be such that $g \in \mathcal{H} \cap \mathcal{C}$ for all $g \in \mathcal{H}$ such that $\|g - f\|_{\mathcal{H}} \le 2\epsilon$ (recall that $f \in \text{int}_{||.||_{\mathcal{H}}}(\mathcal{H} \cap \mathcal{C})$). Now, let $g_F \in \mathcal{H}_F$ with $\|g_F - f|_F\|_{\mathcal{H}_F} \le \epsilon$. By Theorem 6 in [10], there exists $\psi \in \mathcal{H}$ such that $\psi|_F = g_F - f_F$, and $\|\psi\|_{\mathcal{H}} \le \|g_F - f|_F\|_{\mathcal{H}_F} + \epsilon \le 2\epsilon$. This implies that $f + \psi$ belongs to $\mathcal{H} \cap \mathcal{C}$. This in turn implies that $(f + \psi)|_F$ belongs to $\mathcal{H}_F \cap \mathcal{C}_F$, with the same arguments as above. Also, we have $(f+\psi)|_F = f|_F + g_F - f|_F = g_F$, so $g_F \in \mathcal{H}_F \cap \mathcal{C}_F$. Hence, $f|_F \in \text{int}_{||.||_{\mathcal{H}_F}}(\mathcal{H}_F \cap \mathcal{C}_F)$.
With the conditions (H1, F), (H2, F), one can check that the proof of Theorem 1 can be carried out, similarly as in [6]. There are just a few modifications, that we now explain.

The main modification is that in [6], the statement of Lemma 1 and its proof need to be adapted to our context. The adapted lemma is Lemma 6, that we state and prove below.

The second modification is that the reference [6] considers functions indexed on [0, 1] while we consider functions indexed on F. This only entails straightforward changes, since we prove convergence or uniform convergence of functions on F, where the set of knots $S_1^{(m)} \times \dots \times S_d^{(m)}$ is dense in F.

The last modification is that in [6], the set of inequality constraints corresponding to $\mathcal{I}_{F,X,y^{(n)}}$ is of the form

$$ \{ f : [0, 1]^d \to \mathbb{R}; f(x^{(i)}) = y_i, i = 1, \dots, n \}, \quad (46) $$

while the set $\mathcal{I}_{F,X,y^{(n)}}$ in our case is of the form

$$ \left\{ f|_{F} : F \to \mathbb{R}; \sum_{j=1}^{q} \lambda_{i,j} f(a_j) = y_i, i = 1, \dots, n \right\}, \quad (47) $$

where the fixed coefficients $(\lambda_{i,1}, \dots, \lambda_{i,q})_{i=1,\dots,n}$ and the points $a_1, \dots, a_q \in F$ are explicited in the proof of Lemma 6 and come from (18).

The difference between (46) and (47) changes the arguments in the proof of [6] only in the first item after (13) there, the change being straightforward. Thus, from this adaptation of the proof in [6], $\hat{Y}_{S^{(m)}|F}$ converges uniformly on F to the function $Y_{F,\text{opt}}$ defined in the text of Theorem 2.

By Proposition 2, this implies that $\hat{Y}_{S^{(m)}}$ converges uniformly on $[0, 1]^d$ to the function $P_{F \to [0,1]^d}(Y_{F,\text{opt}})$. $\square$

For the next lemmas, notice that the definition of $\pi_{S^{(m)}}(f)$ in (2) for functions in $C([0, 1]^d, \mathbb{R})$ can be extended for functions $f \in C(F, \mathbb{R})$. Indeed, it relies only on the sequence of knots $(t_{\ell}^{(S^{(m)})})_{\ell \in \mathcal{L}_{S^{(m)}}}$, which are included in F.

**Lemma 6.** Consider the setting of Theorem 2. Write

$$ F_m = \{ f \in \mathcal{H}_F : \pi_{S^{(m)}}(f)(x^{(i)}) = y_i, i = 1, \dots, n \}. $$

Let $g \in \mathcal{H}_F \cap \mathcal{I}_{F,X,y^{(n)}}$. Then for $m$ large enough we can define $g_m$ by

$$ g_m = \underset{h \in F_m}{\operatorname{argmin}} \|h - g\|_{\mathcal{H}_F}. $$

Furthermore as $m \to \infty$,

$$ \|g_m - g\|_{\mathcal{H}_F} \to 0. $$

The interpretation of Lemma 6 is that functions in $\mathcal{H}_F$ satisfying the $n$ equality constraints are asymptotically well approximated by their projections on $F_m$. The space $F_m$ is the set of functions in $\mathcal{H}_F$ that satisfy the $n$ equality constraints, when interpolated through $\pi_{S^{(m)}}$. Note that this convergence result applies to functions defined on F and follows from the density of the knots on F.

**Proof of Lemma 6.** We introduce notation for the current subdivision and the current left and right neighbors of $(x_j^{(i)})_{i=1,\dots,n,j=1,\dots,d}$, as in Section (4.1.1). For $i \in \{1, \dots, d\}$ and $m \ge m_0$, let us write $x_{m,j,-}^{(i)} = \max\{u \in S_j^{(m)}; u \le x_j^{(i)}\}$ and $x_{m,j,+}^{(i)} = \min\{u \in S_j^{(m)}; u \ge x_j^{(i)}\}$. Also, if $x_j^{(i)} \notin F_j$, we write $\omega_+^m(x_j^{(i)}) = \frac{x_j^{(i)} - x_{m,j,-}^{(i)}}{x_{m,j,+}^{(i)} - x_{m,j,-}^{(i)}}$ and $\omega_-^m(x_j^{(i)}) = 1 - \omega_+^m(x_j^{(i)})$. If $x_j^{(i)} \in F_i$, we write $\omega_+^m(x_j^{(i)}) = \omega_-^m(x_j^{(i)}) = 1/2$.

We then have, for $i = 1, \dots, n$, for $f \in C(F, \mathbb{R})$,

$$ \pi_{S^{(m)}}(f)(x^{(i)}) = \sum_{\epsilon_1, \dots, \epsilon_d \in \{-, +\}} \left( \prod_{j=1}^{d} \omega_{\epsilon_j}^{m}(x_j^{(i)}) \right) f(x_{m,1,\epsilon_1}^{(i)}, \dots, x_{m,d,\epsilon_d}^{(i)}), \quad (48) $$

from (18) and (19).

Reindexing by a single index the $n2^d$ vertices of hypercubes in (48), let us write $q = n2^d$,

$$ \{a_{m,1}, \dots, a_{m,q}\} \subset F $$
and, for $i = 1, \dots, n,$

$$ \{\lambda_{m,i,1}, \dots, \lambda_{m,i,q}\} \subset \mathbb{R}, $$

such that, for $f \in \mathcal{C}(F, \mathbb{R})$,

$$ \pi_{S^{(m)}}(f)(x^{(i)}) = \sum_{j=1}^{q} \lambda_{m,i,j} f(a_{m,j}). $$

Notice that the $(\lambda_{m,i,j})_{i=1,\dots,n, j=1,\dots,q}$ are the products in (48) (also reindexed) and are thus between 0 and 1. Similarly there exist

$$ \{a_1, \dots, a_q\} \in F $$

and, for $i = 1, \dots, n,$

$$ \{\lambda_{i,1}, \dots, \lambda_{i,q}\} $$

such that, for $f \in \mathcal{C}(F, \mathbb{R})$,

$$ P_{F \to [0,1]^d}(f)(x^{(i)}) = \sum_{j=1}^{q} \lambda_{i,j} f(a_j). \qquad (49) $$

Let us show that, for $m$ large enough, there exist $\gamma_1, \dots, \gamma_n \in \mathbb{R}$ such that

$$ \pi_{S^{(m)}} \left( \sum_{i'=1}^{n} \gamma_{i'} \left( \sum_{j'=1}^{q} \lambda_{m,i',j'} k_F(\cdot, a_{m,j'}) \right) \right) (x^{(i)}) = y_i, \quad i=1, \dots, n. \qquad (50) $$

This is equivalent to

$$ \sum_{j=1}^{q} \lambda_{m,i,j} \left( \sum_{i'=1}^{n} \gamma_{i'} \left( \sum_{j'=1}^{q} \lambda_{m,i',j'} k_F(a_{m,j}, a_{m,j'}) \right) \right) = y_i, \quad i=1, \dots, n. $$

This is equivalent to

$$ \sum_{i'=1}^{n} \gamma_{i'} \sum_{j=1}^{q} \sum_{j'=1}^{q} \lambda_{m,i,j} \lambda_{m,i',j'} k_F(a_{m,j}, a_{m,j'}) = y_i, \quad i=1, \dots, n. $$

With $\gamma = (\gamma_1, \dots, \gamma_n)^{\top}$, this is equivalent to

$$ R_m \gamma = y^{(n)}, $$

where $R_m$ is the $n \times n$ matrix with

$$ (R_m)_{i,i'} = \sum_{j=1}^{q} \sum_{j'=1}^{q} \lambda_{m,i,j} \lambda_{m,i',j'} k_F(a_{m,j}, a_{m,j'}). $$

Consider a GP Z on F with continuous trajectories and covariance function $k_F$. Note that this exists by taking the restriction to F of a GP on $[0, 1]^d$ with continuous trajectories and covariance function k. Then $R_m$ is the covariance matrix of the Gaussian vector

$$ \left( \sum_{j=1}^{q} \lambda_{m,i,j} Z(a_{m,j}) \right)_{i=1,\dots,n} = \left( \pi_{S^{(m)}}(Z)(x^{(i)}) \right)_{i=1,\dots,n}. \qquad (51) $$

Define R as the covariance matrix of the Gaussian vector

$$ \left( \sum_{j=1}^{q} \lambda_{i,j} Z(a_j) \right)_{i=1,\dots,n} = \left( P_{F \to [0,1]^d}(Z)(x^{(i)}) \right)_{i=1,\dots,n}. \qquad (52) $$

We have that, for all the (continuous) trajectories of Z, $\pi_{S^{(m)}}(Z)$ converges uniformly to $P_{F \to [0,1]^d}(Z)$ on $[0, 1]^d$ as $m \to \infty$ from Lemma 7. Hence, the Gaussian vector (51) converges almost surely to the Gaussian vector (52) as $m \to \infty$. Thus, by Gaussianity (see for instance [31, Lemma 1]), the covariance matrix of the Gaussian vector (51) converges as $m \to \infty$ to the covariance matrix of the Gaussian vector (52). Hence $R_m$ goes to R as $m \to \infty$.

Let us show that R is invertible. Up to a re-arrangement in the sum (49), we can assume that $a_1, \dots, a_q$ are two by two distinct (in this case we have removed duplicates and q can be smaller than $n2^d$). From Condition 3, for any $v^{(n)} \in \mathbb{R}^n$, there exists a function g in $E_{S^{(m_0)}}$ such that

$$ (g(x^{(i)}))_{i=1,\dots,n} = v^{(n)}. $$
The function $g$ is equal to $P_{F \to [0,1]^d}(g)$ from Proposition 2, item 2. Hence, with this function $g$ we have from (49)

$$ (P_{F \to [0,1]^d}(g)(x^{(i)}))_{i=1,\dots,n} = \left( \sum_{j=1}^q \lambda_{i,j} g(a_j) \right)_{i=1,\dots,n} = v^{(n)}. $$

This means that the application

$$ (h_1, \dots, h_q) \mapsto \left( \sum_{j=1}^{q} \lambda_{i,j} h_j \right)_{i=1,\dots,n} $$

is surjective. In addition, the covariance matrix of $(Z(a_j))_{j=1}^q$ is invertible because $a_1, \dots, a_q$ are two-by-two distinct. Hence the Gaussian vector

$$ \left( \sum_{j=1}^{q} \lambda_{i,j} Z(a_j) \right)_{i=1,\dots,n} $$

has an invertible covariance matrix. Its covariance matrix is $R$ which shows that $R$ is invertible.

Hence, for $m$ large enough, $R_m$ is invertible and there exist $\gamma_1, \dots, \gamma_n$ such that (50) holds and thus $F_m$ is non-empty and thus $g_m$ is well-defined (from the classical projection theorem). Let us take $m$ large enough such that $R_m$ is invertible for the rest of the proof.

We define the spaces $G_0^m$ and $G_1^m$, respectively, as

$$ G_0^m = \left\{ f \in \mathcal{H}_F : \sum_{j=1}^q \lambda_{m,i,j} f(a_{m,j}) = 0, i = 1, \dots, n \right\} $$

and

$$ G_1^m = \operatorname{span} \left( \sum_{j=1}^{q} \lambda_{m,i,j} k_F(\cdot, a_{m,j}); i = 1, \dots, n \right). $$

For arbitrary $f$ in $F_m$, following the proof of Lemma 1 in [6], we have $F_m = f + G_0^m$ and $g_m = f + P_{G_0^m}(g-f)$, where $P_{G_0^m}$ is the orthogonal projection onto $G_0^m$. Therefore, $g - g_m = g - f - P_{G_0^m}(g-f) \in (G_0^m)^\perp = G_1^m$, where $(G_0^m)^\perp$ is the orthogonal space to $G_1^m$. Then there exist $\beta_1^m, \dots, \beta_n^m \in \mathbb{R}$ such that

$$ g - g_m = \sum_{i'=1}^{n} \beta_{i'}^m \sum_{j'=1}^{q} \lambda_{m,i',j'} k_F(\cdot, a_{m,j'}). \quad (53) $$

Hence for $i = 1, \dots, n$ we have

$$ \left( \sum_{j=1}^{q} \lambda_{m,i,j} g(a_{m,j}) \right) - y_i = \sum_{j=1}^{q} \lambda_{m,i,j} (g - g_m)(a_{m,j}) \\ = \sum_{j=1}^{n} \lambda_{m,i,j} \sum_{i'=1}^{n} \beta_{i'}^m \sum_{j'=1}^{q} \lambda_{m,i',j'} k_F(a_{m,j}, a_{m,j'}). $$

Hence $\beta^m = (\beta_1^m, \dots, \beta_n^m)^\top$ is solution of the system

$$ R_m \beta_m = z_m, $$

where

$$ z_m = \left( \sum_{j=1}^{q} \lambda_{m,i,j} g(a_{m,j}) - y_i \right)_{i=1,\dots,n}. $$

The matrix $R^m$ converges to an invertible matrix as seen before. Furthermore, from Lemma 7,

$$ \sum_{j=1}^{q} \lambda_{m,i,j} g(a_{m,j}) = \pi_{S^{(m)}}(g)(x^{(i)}) \to_{m \to \infty} P_{F \to [0,1]^d}(g)(x^{(i)}) = y_i $$
since $g \in \mathcal{I}_{F,X,y^{(n)}}$. Hence the vector $z_m$ goes to zero. Hence $\beta_1^m, \dots, \beta_n^m$ go to zero as $m \to \infty$. Furthermore, the $\lambda_{m,i,j}$'s are non-negative and bounded by 1 and $q$ is fixed, hence, from (53) and the triangle inequality,

$$ \|g - g_m\|_{\mathcal{H}_F} \le nq \max_{i'=1,\dots,n} \max_{j'=1,\dots,q} |\beta_{i'}^m| |\lambda_{m,i',j'}| \max_{x \in F} \sqrt{k(x, x)} \to_{m \to \infty} 0, $$

which concludes the proof.

**Lemma 7.** Consider the setting of Theorem 2. Let $f \in \mathcal{C}(F, \mathbb{R})$. Then as $m \to \infty$, $\pi_{S^{(m)}}(f) - P_{F \to [0,1]^d}(f) \to 0$, uniformly on $[0, 1]^d$.

**Remark 4.** From Remark 1, we know that

$$ \pi_{S^{(m)}}(f) = P_{F_{S^{(m)}} \to [0,1]^d}(f), $$

where $F_{S^{(m)}} = \prod_{j=1}^d (S_j^{(m)} \cap [0, 1])$. The set $F_{S^{(m)}}$ is contained in $F$ and converges to $F$, for the Hausdorff distance, as $m \to \infty$. Hence, Lemma 7 is a continuity property of the application $F' \mapsto P_{F' \to [0,1]^d}(f)$, for fixed $f$, with respect to $F'$ and for the Hausdorff distance.

**Proof of Lemma 7.** Let $f \in \mathcal{C}(F, \mathbb{R})$. Assume that as $m \to \infty$, $\pi_{S^{(m)}}(f) - P_{F \to [0,1]^d}(f)$ does not go to zero. Then, up to extracting subsequences and by compacity of $[0, 1]^d$, there exist $\epsilon > 0$ and a sequence $(x^{(m)})_{m \ge m_0}$ converging to $x^{(\infty)} \in [0, 1]^d$ such that

$$ \left| \pi_{S^{(m)}}(f)(x^{(m)}) - P_{F \to [0,1]^d}(f)(x^{(m)}) \right| \ge \epsilon. \quad (54) $$

Let us now contradict (54). Similarly to the proof of Lemma 6, we introduce notation for the current subdivision and the current left and right neighbors of $x_1^{(m)}, \dots, x_d^{(m)}$. For $j \in \{1, \dots, d\}$ and $m \ge m_0$, let us write $x_{m,j,-}^{(m)} = \max\{u \in S_j^{(m)}; u \le x_j^{(m)}\}$ and $x_{m,j,+}^{(m)} = \min\{u \in S_j^{(m)}; u \ge x_j^{(m)}\}$. Also, if $x_j^{(m)} \notin F_j$, we write $\omega_+^m(x_j^{(m)}) = \frac{x_j^{(m)} - x_{m,j,-}^{(m)}}{x_{m,j,+}^{(m)} - x_{m,j,-}^{(m)}}$ and $\omega_-^m(x_j^{(m)}) = 1 - \omega_+^m(x_j^{(m)})$. If $x_j^{(m)} \in F_j$, we write $\omega_+^m(x_j^{(m)}) = \omega_-^m(x_j^{(m)}) = 1/2$.

We then have

$$ \pi_{S^{(m)}}(f)(x^{(m)}) = \sum_{\epsilon_1, \dots, \epsilon_d \in \{-, +\}} \left( \prod_{j=1}^{d} \omega_{\epsilon_j}^{m}(x_j^{(m)}) \right) f(x_{m,1,\epsilon_1}^{(m)}, \dots, x_{m,d,\epsilon_d}^{(m)}), $$

from (18) and (19). Up to extracting subsequences, we can partition $\{1, \dots, d\}$ as $J_1 \cup J_2 \cup J_3 \cup J_4$ (where some of the four sets are possibly empty) where

* for $j \in J_1$ we have $x_j^{(\infty)} \notin F_j$,

* for $j \in J_2$ we have $x_j^{(\infty)} \in F_j$, $x_{m,j,-}^{(m)} \to x_j^{(\infty)}$ and $x_{m,j,+}^{(m)} \to x_j^{(\infty)}$,

* for $j \in J_3$ we have $x_j^{(\infty)} \in F_j$, $x_{m,j,-}^{(m)} \to x_j^{(\infty)}$ and $\liminf |x_{m,j,+}^{(m)} - x_j^{(\infty)}| > 0$,

* and for $j \in J_4$ we have $x_j^{(\infty)} \in F_j$, $\liminf |x_{m,j,-}^{(m)} - x_j^{(\infty)}| > 0$ and $x_{m,j,+}^{(m)} \to x_j^{(\infty)}$.

Note that we can indeed find this partition, after extraction of subsequences, because for $x_j^{(\infty)} \in F_j$, we have $\liminf |x_{m,j,-}^{(m)} - x_j^{(\infty)}| = 0$ or $\liminf |x_{m,j,+}^{(m)} - x_j^{(\infty)}| = 0$. Up to re-indexing and without loss of generality, we assume that $J_1 = \{1, \dots, j_1\}$, $J_2 = \{j_1 + 1, \dots, j_2\}$, $J_3 = \{j_2 + 1, \dots, j_3\}$ and $J_4 = \{j_3 + 1, \dots, d\}$.

We have

$$ \left| \pi_{S^{(m)}}(f)(x^{(m)}) - P_{F \to [0,1]^d}(f)(x^{(m)}) \right| = \left| \sum_{\epsilon_1, \dots, \epsilon_d \in \{-, +\}} \left( \prod_{j=1}^{d} \omega_{\epsilon_j}^{m}(x_j^{(m)}) \right) f(x_{m,1,\epsilon_1}^{(m)}, \dots, x_{m,d,\epsilon_d}^{(m)}) - \sum_{\epsilon_1, \dots, \epsilon_d \in \{-, +\}} \left( \prod_{j=1}^{d} \omega_{\epsilon_j}(x_j^{(m)}) \right) f(x_{1,\epsilon_1}^{(m)}, \dots, x_{d,\epsilon_d}^{(m)}) \right|. \quad (55) $$
In the above display, for $j \in J_3$, then $\liminf |x_{m,j,+}^{(m)} - x_j^{(\infty)}| > 0$ and $x_{m,j,-}^{(m)} \to x_j^{(\infty)}$. Hence one can see that if $\epsilon_j = +$, $w_{\epsilon_j}^m(x_j^{(m)}) \to 0$. Similarly, for $j \in J_4$, if $\epsilon_j = -$, then $w_{\epsilon_j}^m(x_j^{(m)}) \to 0$.

Similarly, for $j \in J_3$, if $\epsilon_j = +$, then $w_{\epsilon_j}^m(x_j^{(m)}) \to 0$ and for $j \in J_4$, if $\epsilon_j = -$, then $w_{\epsilon_j}^m(x_j^{(m)}) \to 0$.

Hence, in the two sums in (55), if $(\epsilon_1, \dots, \epsilon_d)$ is such that $\epsilon_j = +$ for $j \in J_3$, or $\epsilon_j = -$ for $j \in J_4$, then the corresponding summand goes to zero. As a consequence we have

$$ \left| \pi_{S^{(m)}}(f)(x^{(m)}) - P_{F \to [0,1]^d}(f)(x^{(m)}) \right| = o(1) + \left| \sum_{\epsilon_1, \dots, \epsilon_{j_2} \in \{-, +\}} \left( \prod_{j=1}^{j_2} \omega_{\epsilon_j}^m(x_j^{(m)}) \right) f(x_{m,1,\epsilon_1}^{(m)}, \dots, x_{m,j_2,\epsilon_{j_2}}^{(m)}, \right. \\ \left. \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad x_{m,j_2+1,-}^{(m)}, \dots, x_{m,j_3,-}^{(m)}, x_{m,j_3+1,+}^{(m)}, \dots, x_{m,j_4,+}^{(m)} \right. \\ \left. \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad - \sum_{\epsilon_1, \dots, \epsilon_{j_2} \in \{-, +\}} \left( \prod_{j=1}^{j_2} \omega_{\epsilon_j}(x_j^{(m)}) \right) f(x_{1,\epsilon_1}^{(m)}, \dots, x_{j_2,\epsilon_{j_2}}^{(m)}, x_{j_2+1,-}^{(m)}, \dots, x_{j_3,-}^{(m)}, x_{j_3+1,+}^{(m)}, \dots, x_{j_4,+}^{(m)}) \right| $$

For $j \in \{j_1 + 1, \dots, j_2\}$, we have $x_j^{(\infty)} \in F_j$, $x_{m,j,-}^{(m)} \to x_j^{(\infty)}$ and $x_{m,j,+}^{(m)} \to x_j^{(\infty)}$. As a consequence, because $x_{m,j,-}^{(m)} \le x_{j,-}^{(m)} \le x_j^{(m)} \le x_{j,+}^{(m)} \le x_{m,j,+}^{(m)}$, also $x_{j,-}^{(m)} \to x_j^{(\infty)}$ and $x_{j,+}^{(m)} \to x_j^{(\infty)}$. Similarly, for $j \in \{j_2 + 1, \dots, j_3\}$, we have $x_{m,j,-}^{(m)} \to x_j^{(\infty)}$ and so also $x_{j,-}^{(m)} \to x_j^{(\infty)}$. Similarly, for $j \in \{j_3 + 1, \dots, d\}$, we have $x_{m,j,+}^{(m)} \to x_j^{(\infty)}$ and so also $x_{j,+}^{(m)} \to x_j^{(\infty)}$.

This yields, by continuity of $f$,

$$ \left| \pi_{S^{(m)}}(f)(x^{(m)}) - P_{F \to [0,1]^d}(f)(x^{(m)}) \right| = o(1) + \left| \sum_{\epsilon_1, \dots, \epsilon_{j_2} \in \{-, +\}} \left( \prod_{j=1}^{j_2} \omega_{\epsilon_j}^m(x_j^{(m)}) \right) f(x_{m,1,\epsilon_1}^{(m)}, \dots, x_{m,j_1,\epsilon_{j_1}}^{(m)}, x_{j_1+1,\epsilon_1}^{(\infty)}, \dots, x_{j_2}^{(\infty)}, x_{j_2+1,\epsilon_1}^{(\infty)}, \dots, x_{j_3}^{(\infty)}, x_{j_3+1,\epsilon_1}^{(\infty)}, \dots, x_{j_4}^{(\infty)}) \right. \\ \left. - \sum_{\epsilon_1, \dots, \epsilon_{j_2} \in \{-, +\}} \left( \prod_{j=1}^{j_2} \omega_{\epsilon_j}(x_j^{(m)}) \right) f(x_{1,\epsilon_1}^{(m)}, \dots, x_{j_1,\epsilon_{j_1}}^{(m)}, x_{j_1+1,\epsilon_1}^{(\infty)}, \dots, x_{j_2}^{(\infty)}, x_{j_2+1,\epsilon_1}^{(\infty)}, \dots, x_{j_3}^{(\infty)}, x_{j_3+1,\epsilon_1}^{(\infty)}, \dots, x_{j_4}^{(\infty)}) \right| $$

In the above display, for the first sum, we can isolate

$$ \sum_{\epsilon_{j_1+1}, \dots, \epsilon_{j_2} \in \{-, +\}} \left( \prod_{j=j_1+1}^{j_2} \omega_{\epsilon_j}^m(x_j^{(m)}) \right) = 1, $$

because the argument of $f$ does not depend on $\epsilon_{j_1+1}, \dots, \epsilon_{j_2}$. We can proceed similarly with the second sum. This yields

$$ \left| \pi_{S^{(m)}}(f)(x^{(m)}) - P_{F \to [0,1]^d}(f)(x^{(m)}) \right| = o(1) + \left| \sum_{\epsilon_1, \dots, \epsilon_{j_1} \in \{-, +\}} \left( \prod_{j=1}^{j_1} \omega_{\epsilon_j}^m(x_j^{(m)}) \right) f(x_{m,1,\epsilon_1}^{(m)}, \dots, x_{m,j_1,\epsilon_{j_1}}^{(m)}, x_{j_1+1,\epsilon_1}^{(\infty)}, \dots, x_{j_2}^{(\infty)}, x_{j_2+1,\epsilon_1}^{(\infty)}, \dots, x_{j_3}^{(\infty)}, x_{j_3+1,\epsilon_1}^{(\infty)}, \dots, x_{j_4}^{(\infty)}) \right. \\ \left. - \sum_{\epsilon_1, \dots, \epsilon_{j_1} \in \{-, +\}} \left( \prod_{j=1}^{j_1} \omega_{\epsilon_j}(x_j^{(m)}) \right) f(x_{1,\epsilon_1}^{(m)}, \dots, x_{j_1,\epsilon_{j_1}}^{(m)}, x_{j_1+1,\epsilon_1}^{(\infty)}, \dots, x_{j_2}^{(\infty)}, x_{j_2+1,\epsilon_1}^{(\infty)}, \dots, x_{j_3}^{(\infty)}, x_{j_3+1,\epsilon_1}^{(\infty)}, \dots, x_{j_4}^{(\infty)}) \right| $$
Now for $j = 1, \dots, j_1$, we have $x_{j,-}^{(\infty)} < x_j^{(\infty)} < x_{j,+}^{(\infty)}$. Hence by continuity, we have

$$
\begin{aligned}
& \left| \pi_{S^{(m)}}(f)(x^{(m)}) - P_{F \to [0,1]^d}(f)(x^{(m)}) \right| \\
& = o(1) + \left| \sum_{\epsilon_1, \dots, \epsilon_{j_1} \in \{-, +\}} \left( \prod_{j=1}^{j_1} \omega_{\epsilon_j}(x_j^{(\infty)}) \right) \right. \\
& \qquad \left. f(x_{1,\epsilon_1}^{(\infty)}, \dots, x_{j_1,\epsilon_{j_1}}^{(\infty)}, x_{j_1+1}^{(\infty)}, \dots, x_{j_2}^{(\infty)}, x_{j_2+1}^{(\infty)}, \dots, x_{j_3}^{(\infty)}, x_{j_3+1}^{(\infty)}, \dots, x_{j_4}^{(\infty)}) \right. \\
& \qquad \left. - \sum_{\epsilon_1, \dots, \epsilon_{j_1} \in \{-, +\}} \left( \prod_{j=1}^{j_1} \omega_{\epsilon_j}(x_j^{(\infty)}) \right) \right. \\
& \qquad \left. f(x_{1,\epsilon_1}^{(\infty)}, \dots, x_{j_1,\epsilon_{j_1}}^{(\infty)}, x_{j_1+1}^{(\infty)}, \dots, x_{j_2}^{(\infty)}, x_{j_2+1}^{(\infty)}, \dots, x_{j_3}^{(\infty)}, x_{j_3+1}^{(\infty)}, \dots, x_{j_4}^{(\infty)}) \right|.
\end{aligned}
$$

Notice that the two sums above are the same. Hence

$$
\left| \pi_{S^{(m)}}(f)(x^{(m)}) - P_{F \to [0,1]^d}(f)(x^{(m)}) \right| = o(1).
$$

This is in contradiction with (54) which concludes the proof.

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

**Proof of Theorem 3.** Since the sequence of sets $\mathcal{J}_m$ is nested, it is equal to a set $\mathcal{J}_\infty$ for $m$ larger than some $m_0 \in \mathbb{N}$. Let us consider $m \ge m_0$ for the rest of the proof. Let $d = |\mathcal{J}_\infty|$. Then, with the set of variables $\mathcal{J}_\infty$, and the sequence of subdivisions $(S^{(m)})_{m \ge m_0}$ in $S_{\mathcal{J}_\infty}$, we can check that the conditions of Theorem 3 imply the conditions of Theorem 2.

Hence, from Theorem 2, as $m \to \infty$, the sequence of functions $\hat{Y}_{\mathcal{J}_m, S^{(m)}}$ converges uniformly to a limit function. Hence, as $m \to \infty$, we have,

$$
I_{\mathcal{J}_m, S^{(m)}}(i_{m+1}^*, t_{m+1}^*) \to 0.
$$

Assume now that there exists $i \in \mathcal{J}_\infty$ such that $S_i^{(m)}$ is not dense in $[0, 1]$. As for $t_i \in [0, 1]$, the sequence $d(t_i, S_i^{(m)})$ is decreasing (with respect to $m$) and thus has a limit, this implies that there exists $\epsilon > 0$ and $t_i \in [0, 1]$ such that $d(t_i, S_i^{(m)}) \ge \epsilon$ for all $m \in \mathbb{N}$. This means that for $m$ large enough such that $b_m \le \epsilon$,

$$
\sup_{\substack{i \in \mathcal{J}_m, t \in [0,1], \\ d(t, S_i^{(m)}) \ge b_m}} \left( I_{\mathcal{J}_m, S^{(m)}}(i,t) + \Delta d(t, S_i^{(m)}) \right) \ge \Delta \epsilon.
$$

As $m \ge m_0$, then $i_{m+1}^* \in \mathcal{J}_m = \mathcal{J}_\infty$. Thus, by definition of algorithm 1,

$$
I_{\mathcal{J}_m, S^{(m)}}(i_{m+1}^*, t_{m+1}^*) + \Delta d(t_{m+1}^*, S_{i_{m+1}^*}^{(m)}) + a_m \geq \Delta \epsilon.
$$

Hence, we will reach a contradiction if we show that $\liminf_{m \to \infty} d(t_{m+1}^*, S_{i_{m+1}^*}^{(m)}) = 0$. There is at least one coordinate $i \in \mathcal{J}_m$ chosen an infinite number of times by MaxMod algorithm. Let $(m_\ell)_{\ell \in \mathbb{N}}$ be the corresponding subsequence of values of $m$, i.e. for which $i_{m+1}^* = i$. Then we have $\{t_{m_1+1}^*, \dots, t_{m_\ell+1}^*\} \subset S_i^{(m_\ell+1)}$ and thus $d(t_{m_1+1}^*, S_i^{(m_\ell+1)}) \le d(t_{m_\ell+1}^*, \{t_{m_1}^*, \dots, t_{m_\ell}^*\})$. The limit inferior of this last quantity is zero, by considering a convergent subsequence of $(t_{m_\ell}^*)_{\ell \in \text{the compact interval } [0, 1]}$. Hence, we have reached a contradiction and, eventually, for all $i \in \{1, \dots, d\}$, $S_i^{(m)}$ is dense in $[0, 1]$.

Let us now assume that $d < D$. Let us then consider $j_\infty \in \{1, \dots, D\} \setminus \mathcal{J}$. We then have by definition of the MaxMod algorithm, for $m$ larger than $m_0$,

$$
I_{\mathcal{J}_m, S^{(m)}}(i_{m+1}^*, t_{m+1}^*) + \Delta d(t_{m+1}^*, S_{i_{m+1}^*}^{(m)}) + a_m \geq \Delta',
$$

by considering $i = j_\infty \notin \mathcal{J}_\infty$ in (13). The left hand side of the above display goes to zero, because for $i \in \mathcal{J}_\infty$, $S_i^{(m)}$ is dense in $[0, 1]$, as we have shown and also because

$$
I_{\mathcal{J}_m, S^{(m)}}(i_{m+1}^*, t_{m+1}^*) \to 0
$$

as we have shown. This yields a contradiction. Hence $\mathcal{J}_m = \{1, \dots, D\}$ for $m \ge m_0$.

From the fact that all the variables are eventually active and from the density of $S_i^{(m)}$ for $i \in \{1, \dots, D\}$, we conclude from Theorem 1 (Theorem 3.2 in [6]). $\square$
References

[1] D. Azzimonti. profExtrema: Compute and visualize profile extrema functions. https://cran.r-project.org/web/packages/profExtrema/index.html, 2018.

[2] D. Azzimonti, D. Ginsbourger, J. Rohmer, and D. Idier. Profile extrema for visualizing and quantifying uncertainties on excursion regions. Application to coastal flooding. Technometrics, 61(4):474–493, 2019.

[3] F. Bachoc, A. Lagnoux, and A. F. López-Lopera. Maximum likelihood estimation for Gaussian processes under inequality constraints. Electronic Journal of Statistics, 13(2):2921–2969, 2019.

[4] X. Bay and B. Gauthier. Extended formula for Kriging interpolation. Technical report, Ecole des Mines de St-Etienne, 2010.

[5] X. Bay, L. Grammont, and H. Maatouk. Generalization of the Kimeldorf-Wahba correspondence for constrained interpolation. Electronic journal of statistics, 10(1):1580–1595, 2016.

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

[7] J. Bect, F. Bachoc, and D. Ginsbourger. A supermartingale approach to Gaussian process based sequential design of experiments. Bernoulli, 25(4A):2883–2919, 2019.

[8] P. Bellec. Sharp oracle inequalities for least squares estimators in shape restricted regression. The Annals of Statistics, 46(2):745–780, 2018.

[9] M. Ben Salem, O. Roustant, F. Gamboa, and L. Tomaso. Universal prediction distribution for surrogate models. SIAM/ASA Journal on Uncertainty Quantification, 5(1):1086–1109, 2017.

[10] A. Berlinet and C. Thomas-Agnan. Reproducing kernel Hilbert spaces in probability and statistics. Springer Science & Business Media, 2011.

[11] S. Boyd, S. P. Boyd, and L. Vandenberghe. Convex optimization. Cambridge university press, 2004.

[12] J.-P. Chiles and P. Delfiner. Geostatistics: Modeling Spatial Uncertainty. John Wiley & Sons, 2009.

[13] A. Cousin, H. Maatouk, and D. Rullière. Kriging of financial term-structures. European Journal of Operational Research, 255(2):631–648, 2016.

[14] J. Creutzig, T. Müller-Gronbach, and K. Ritter. Free-knot spline approximation of stochastic processes. Journal of Complexity, 23(4):867–889, 2007.

[15] S. Da Veiga and A. Marrel. Gaussian process modeling with inequality constraints. Annales de la faculté des sciences de Toulouse Mathématiques, 21(3):529–555, 2012.

[16] S. Da Veiga and A. Marrel. Gaussian process regression with linear inequality constraints. Reliability Engineering & System Safety, 195:106732, 2020.

[17] C. De Boor. A Practical Guide to Splines. Applied Mathematical Sciences. Springer, New York, 2001.

[18] C. De Boor. Spline basics. In Handbook of Computer Aided Geometric Design, pages 141–163. North-Holland, Amsterdam, 2002.

[19] D. Dupuy, C. Helbert, and J. Franco. DiceDesign and DiceEval: Two R packages for design and analysis of computer experiments. Journal of Statistical Software, 65(11):1–38, 2015.

[20] C. Durot. Sharp asymptotics for isotonic regression. Probability theory and related fields, 122(2):222–240, 2002.

[21] C. Durot and H. P. Lopuhaä. Limit theory in monotone function estimation. Statistical Science, 33(4):547–567, 2018.
[22] H. G. Eggleston. *Convexity*. Cambridge Tracts in Mathematics. Cambridge University Press, 1958.

[23] M. G. Genton. Classes of kernels for machine learning: A statistics perspective. *Journal of Machine Learning Research*, 2:299–312, 2001.

[24] S. Golchi, D. R. Bingham, H. Chipman, and D. A. Campbell. Monotone emulation of computer experiments. *SIAM/ASA Journal on Uncertainty Quantification*, 3(1):370–392, 2015.

[25] D. Goldfarb and A. Idnani. A numerically stable dual method for solving strictly convex quadratic programs. *Mathematical programming*, 27(1):1–33, 1983.

[26] R. Goldman. B-spline approximation and the de Boor algorithm. In *Pyramid Algorithms*, pages 347–443. Elsevier, Amsterdam, 2003.

[27] P. Groeneboom and G. Jongbloed. *Nonparametric estimation under shape constraints*, volume 38. Cambridge University Press, 2014.

[28] P. Groeneboom, G. Jongbloed, and J. Wellner. Estimation of a convex function: characterizations and asymptotic theory. *Annals of Statistics*, 29(6):1653–1698, 2001.

[29] U. Hornung. Monotone spline-interpolation. In *Numerische Methoden der Approximationstheorie*, pages 172–191. Springer, 1978.

[30] Y. Hu. An algorithm for data reduction using splines with free knots. *IMA Journal of Numerical Analysis*, 13(3):365–381, 1993.

[31] I. A. Ibragimov and Y. A. Rozanov. *Gaussian random processes*. Springer Science & Business Media, 2012.

[32] D. R. Jones, M. Schonlau, and W. J. Welch. Efficient global optimization of expensive black-box functions. *Journal of Global Optimization*, 13(4):455–492, Dec 1998.

[33] D. Jupp. Approximation to data by splines with free knots. *SIAM Journal on Numerical Analysis*, 15(2):328–343, 1978.

[34] M. C. Kennedy and A. O'Hagan. Bayesian calibration of computer models. *Journal of the Royal Statistical Society: Series B (Statistical Methodology)*, 63(3):425–464, 2001.

[35] L. P. Kobbelt. Multiresolution techniques. In *Handbook of Computer Aided Geometric Design*, pages 343–361. North-Holland, Amsterdam, 2002.

[36] L. Lin and D. B. Dunson. Bayesian monotone regression using Gaussian process projection. *Biometrika*, 101(2):303–317, 2014.

[37] A. F. López-Lopera. *lineqGPR: Gaussian Process Regression Models with Linear Inequality Constraints*, 2019. R package version 0.2.0.

[38] A. F. López-Lopera, F. Bachoc, N. Durrande, J. Rohmer, D. Idier, and O. Roustant. Approximating Gaussian process emulators with linear inequality constraints and noisy observations via MC and MCMC. In *Monte Carlo and Quasi-Monte Carlo Methods*, pages 363–381. Springer International Publishing, 2020.

[39] A. F. López-Lopera, F. Bachoc, N. Durrande, and O. Roustant. Finite-dimensional Gaussian approximation with linear inequality constraints. *SIAM/ASA Journal on Uncertainty Quantification*, 6(3):1224–1255, 2018.

[40] A. F. López-Lopera, S. John, and N. Durrande. Gaussian process modulated Cox processes under linear inequality constraints. In *International Conference on Artificial Intelligence and Statistics*, pages 1997–2006, 2019.

[41] H. Maatouk. Finite-dimensional approximation of Gaussian processes with inequality constraints and errors measurements. *Preprint HAL*, 2017.

[42] H. Maatouk and X. Bay. Gaussian process emulators for computer experiments with inequality constraints. *Mathematical Geosciences*, 49(5):557–582, 2017.
[43] M. Niu, P. Cheung, L. Lin, Z. Dai, N. Lawrence, and D. Dunson. Intrinsic Gaussian processes on complex constrained domains. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 81(3):603–627, 2019.

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

[45] C. E. Rasmussen and C. K. I. Williams. Gaussian Processes for Machine Learning. The MIT Press, Cambridge, MA, 2005.

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

[47] J. Riihimäki and A. Vehtari. Gaussian processes with monotonicity information. In International Conference on Artificial Intelligence and Statistics, pages 645–652, 2010.

[48] J. Rohmer, D. Idier, F. Paris, R. Pedreros, and J. Louisor. 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, 2018.

[49] J. Sacks, W. Welch, T. Mitchell, and H. Wynn. Design and analysis of computer experiments. Statistical Science, 4:409–423, 1989.

[50] M. Slassi. The optimal free knot spline approximation of stochastic differential equations with additive noise. Journal of Computational and Applied Mathematics, 261:62–71, 2014.

[51] M. Stein. Interpolation of Spatial Data: Some Theory for Kriging. Springer Series in Statistics. Springer, New York, 1999.

[52] S. Zhou, P. Giulani, J. Piekarewicz, A. Bhattacharya, and D. Pati. Reexamining the proton-radius problem using constrained Gaussian processes. Physical Review C, 99:055202, 2019.