Nested Kriging predictions for datasets with a large number of observations

This work falls within the context of predicting the value of a real function at some input locations given a limited number of observations of this function. The Kriging interpolation technique (or Gaussian process regression) is often considered to tackle such a problem, but the method suffers from its computational burden when the number of observation points is large. We introduce in this article nested Kriging predictors which are constructed by aggregating sub-models based on subsets of observation points. This approach is proven to have better theoretical properties than other aggregation methods that can be found in the literature. Contrarily to some other methods it can be shown that the proposed aggregation method is consistent. Finally, the practical interest of the proposed method is illustrated on simulated datasets and on an industrial test case with \(10^4\) observations in a 6-dimensional space.

Part of this research was conducted within the frame of the Chair in Applied Mathematics OQUAIDO, gathering partners in technological research (BRGM, CEA, IFPEN, IRSN, Safran, Storengy) and academia (Ecole Centrale de Lyon, Mines Saint-Etienne, University of Grenoble, University of Nice, University of Toulouse and CNRS) around advanced methods for Computer Experiments. The authors would like to warmly thank Dr. Géraud Blatman and EDF R&D for providing us the industrial test case. They also thank both editor and reviewers for very precise and constructive comments on this paper. This paper has been finished during a stay of D. Rullière at Vietnam Institute for Advanced Study in Mathematics, the latter author thanks the VIASM institute and DAMI research chair (Data Analytics & Models for Insurance) for their support.

Authors and Affiliations

1. Laboratoire SAF, EA2429, ISFA, Université Claude Bernard Lyon 1, Université de Lyon, 50 avenue Tony Garnier, 69366, Lyon, France

   Didier Rullière

2. Institut Fayol—LIMOS, Mines Saint-Étienne, 158 cours Fauriel, Saint-Étienne, France

   Nicolas Durrande

3. Institut de Mathématiques de Toulouse, Université Paul Sabatier, 118 route de Narbonne, 31062, Toulouse, France

   François Bachoc

4. Institute of Statistics, University of Neuchâtel, avenue de Bellevaux 51, 2000, Neuchâtel, Switzerland

   Clément Chevalier

Corresponding author

Correspondence to Nicolas Durrande.

Appendix: Proof of Proposition 4

Complexities Under chosen assumption on \(\alpha \) and \(\beta \) coefficients, for a regular tree and in the case of simple Kriging sub-models, \(\mathcal {C}_\alpha =\sum _{\nu =1}^{\bar{\nu }} \sum _{i=1}^{n_\nu } \alpha c_\nu ^3 =\alpha \sum _{\nu =1}^{\bar{\nu }} c_\nu ^3 n_\nu \) and \(\mathcal {C}_\beta =\sum _{\nu =1}^{\bar{\nu }} \sum _{i=2}^{n_\nu } \sum _{j=1}^{i-1}\beta c^2_{\nu } =\frac{\beta }{2} \sum _{\nu =1}^{\bar{\nu }} n_\nu (n_{\nu }-1) c^2_{\nu }\). Notice that the sum starts from \(\nu =1\) in order to include sub-models calculation. Equilibrated trees complexities in a constant child number setting, when \(c_\nu =c\) for all \(\nu \), the tree structure ensures that \(n_{\nu }=n/c^{\nu }\), thus as \(c=n^{1/\bar{\nu }}\), we get when \(n \rightarrow +\infty \), \(\mathcal {C}_\alpha \sim \alpha n^{1+\frac{2}{\bar{\nu }}}\) and \(\mathcal {C}_\beta \sim \frac{\beta }{2} n^2\). The result for equilibrated two-layer tree where \(\bar{\nu }=2\) directly derives from this one, and in this case \(\mathcal {C}_\alpha \sim \alpha n^{2}\) and \(\mathcal {C}_\beta \sim \frac{\beta }{2} n^2\) (it derives also from the expressions of \(\mathcal {C}_\alpha \), \(\mathcal {C}_\beta \), when \(c_1=c_2=\sqrt{n}\), \(n_1=\sqrt{n}\), \(n_2=1\)). Optimal tree complexities one easily shows that under the chosen assumptions \(\mathcal {C}_\beta \sim \frac{\beta }{2}n^2\). Thus, it is indeed not possible to reduce the whole complexity to orders lower than \(O(n^2)\). However, one can choose the tree structure in order to reduce the complexity \(\mathcal {C}_\alpha \). For a regular tree, \(n_\nu =n/(c_1 \cdots c_{\nu })\) such that \(\frac{\partial }{\partial c_k} n_{\nu } = -\mathbf {1}_{{\lbrace {\nu \ge k}\rbrace }} n_\nu /c_k\). Using a Lagrange multiplier \(\ell \), one defines \(\xi (k)=c_k \frac{\partial }{\partial c_k} \left( \mathcal {C}_\alpha - \ell (c_1 \cdots c_{\bar{\nu }} -n) \right) = 3\alpha c_k^3 n_k - \alpha \sum _{\nu =k}^{\bar{\nu }}c_\nu ^3 n_\nu - \ell c_1 \cdots c_{\bar{\nu }}\). The tree structure that minimizes \(\mathcal {C}_\alpha \) is such that for all \(k<\bar{\nu }\), \(\xi (k)=\xi (k+1)=0\). Using \(c_{k+1} n_{k+1}=n_k\), one gets \(3c_{k+1}^2=2 c_{k}^3\) for all \(k<\bar{\nu }\), and setting \(c_1\cdots c_{\bar{\nu }}=n\), \(c_\nu = \delta \left( \delta ^{-\bar{\nu }} n\right) ^{\frac{\delta ^{\nu -1}}{2(\delta ^{\bar{\nu }}-1)}}\), \(\nu =1, \ldots , \bar{\nu }\), with \(\delta =\frac{3}{2}\). Setting \(\gamma =\frac{27}{4}\delta ^{-\frac{\bar{\nu }}{\delta ^{\bar{\nu }}-1}}\left( 1-\delta ^{-\bar{\nu }}\right) \). After some direct calculations this tree structure corresponds to complexities, \(\mathcal {C}_\alpha = \gamma \alpha n^{1+\frac{1}{\delta ^{\bar{\nu }}-1}}\) and \(\mathcal {C}_\beta \sim \frac{\beta }{2}n^2\). In a two-layers setting one gets \(c_1=\left( \frac{3}{2}\right) ^{1/5} n^{2/5}\) and \(c_2=\left( \frac{3}{2}\right) ^{-1/5} n^{3/5}\), which leads to \(\mathcal {C}_\alpha = \gamma \alpha n^{9/5}\) and \(\mathcal {C}_\beta = \frac{\beta }{2} n^2 - \frac{\beta }{2} \left( \frac{3}{2}\right) ^{\frac{1}{5}}n^{\frac{7}{5}}\), where \(\gamma =(\frac{2}{3})^{-2/5}+(\frac{2}{3})^{3/5}\simeq 1.96\) (eventually notice that even for values of n of order \(10^5\), terms of order like \(n^{9/5}\) are not necessarily negligible compared to those of order \(n^2\), and that \(\mathcal {C}_\beta \) is slightly affected by the choice of the tree structure, but the global complexity benefits from the optimization of \(\mathcal {C}_\alpha \)).

Storage footprint First, covariances can be stored in triangular matrices. So temporary objects M, k and K in Algorithm 1 require the storage of \(c_{\max }(c_{\max }+5)/2\) real values. For a given step \(\nu \), \(\nu \ge 2\), building all vectors \(\alpha _i\) requires the storage of \(\sum _{i=1}^{n_{\nu }} c_i^{\nu }=n_{\nu -1}\) values. At last, for a given step \(\nu \), we simultaneously need objects \(M_{\nu -1}, K_{\nu -1}, M_{\nu }, K_{\nu }\), which require the storage of \(n_{\nu -1}(n_{\nu -1}+3)/2 + n_{\nu }(n_{\nu }+3)/2\) real values. In a regular tree, as \(n_\nu \) is decreasing in \(\nu \), the storage footprint is \(\mathcal {S} = (c_{\max }(c_{\max }+5) + n_1(n_1+5) + n_2(n_2+3))/2\). Hence the equivalents for \(\mathcal {S}\) for the different tree structures, \(\mathcal {S}\sim n\) for the two-layer equilibrated tree, \(\mathcal {S}\sim \frac{1}{2}n^{2-2/\bar{\nu }}\) for the \(\bar{\nu }\)-layer, \(\bar{\nu }>2\) and the indicated result for the optimal tree. Simple orders are given in the proposition, which avoids separating the case \(\bar{\nu }=2\) and a cumbersome constant for the optimal tree.

Reprints and permissions