A criterion and incremental design construction for simultaneous kriging predictions

In this paper, we further investigate the problem of selecting a set of design points for universal kriging, which is a widely used technique for spatial data analysis. Our goal is to select the design points in order to make simultaneous predictions of the random variable of interest at a finite number of unsampled locations with maximum precision. Specifically, we consider as response a correlated random field given by a linear model with an unknown parameter vector and a spatial error correlation structure. We propose a new design criterion that aims at simultaneously minimizing the variation of the prediction errors at various points. We also present various efficient techniques for incrementally building designs for that criterion scaling well for high dimensions. Thus the method is particularly suitable for big data applications in areas of spatial data analysis such as mining, hydrogeology, natural resource monitoring, and environmental sciences or equivalently for any computer simulation experiments. We have demonstrated the effectiveness of the proposed designs through two illustrative examples: one by simulation and another based on real data from Upper Austria.


Introduction
Given a finite number k +m of locations x we are interested into making simultaneous predictions Ŷ (•) of Y (•) at m unsampled locations using observations Y (x 1 ), . . ., Y (x k ) collected at some design points ξ = (x 1 , . . ., x k ) ⊂ X k .Our objective is to select ξ (of given size k) in order to maximize the precision of the predictions Ŷ (x) over X .This setup is used in such diverse areas of spatial data analysis as mining, hydrogeology, natural resource monitoring and environmental sciences, see, e.g., Cressie (1993), and has become the standard modeling paradigm in computer simulation experiments (cf.Fang et al. (2005); Kleijnen (2009); Rasmussen and Williams (2005); Santner et al. (2003)), known under the designations of Gaussian Process (GP) modelling and kriging analysis.For a general review in the context of spatial statistics see Wang et al. (2012).
Specifically, the model underlying our investigations is the model for universal kriging, i.e. we have a correlated scalar random field given by Here, β is an unknown vector of parameters in R p , µ(•, •) a known function of regressors at some given locations x in a compact subset X of R d and the random term ε (x) has zero mean, variance σ 2 and a parameterized spatial error correlation structure such that IE (ε (x) ε (x ′ )) = σ 2 c(x, x ′ ; ν) with ν some covariance parameters.We further assume that • the deterministic term is linear in the parameters β, i.e., µ(x, β) = f(x)β, where ) is a vector of known functions, • the first two moments of the error ε (x) and hence of Y (x) exist and • the variance σ 2 and the covariance parameters ν are known.
It is often assumed that the random field ε (x) is Gaussian, allowing estimation of β and θ = {σ 2 , ν} by Maximum Likelihood.We do not need to assume a stationary nor an isotropic covariance structure.
Traditional optimality criteria for designs for prediction are kriging variance-minimizing, albeit not all authors understand the same by the term kriging variance or kriging covariance.Some mean the variance of the prediction Var Ŷ (x) , cf. (Müller et al. (2015)), some the variance of the prediction error Var Ŷ (x) − Y (x) , cf. (Cressie (1993)), and some even Var (Y (x)|Y (ξ)), i.e. the variance of Y (•) at unsampled locations x given the observations on the design points ξ, cf.(Chevalier and Ginsbourger (2012)).We follow the second perception because trying to minimize the variation of the prediction errors seems to yield most precise predictions.Definition 1.Let x ∈ X be an arbitrary unsampled location and Ŷ (x) the best linear unbiased predictor (BLUP) at x.The kriging variance at x is σ 2 (x) := Var Ŷ (x) − Y (x) i.e. the variance of the best linear predictor minus the random variable to be predicted.
Let x ′ ̸ = x be a second unsampled location and Ŷ (x ′ ) the BLUP at x ′ ∈ X .The kriging covariance for x and x ′ is With the above definition G-optimal designs (cf.Dasgupta et al. (2022a)) try to minimize the maximum kriging variance, i.e. (2) Another popular optimality criterion tries to minimize the average prediction variance over a set of m specific points (V-optimality), i.e.
Note that the latter if not supported on a grid but rather covering the whole X , expressed as an integral is often called I-optimality, see Dasgupta et al. (2022b) for a recent example of the terminology.
None of these and other criteria for designs for prediction considers the kriging covariances or the kriging covariance matrix at the unsampled locations (x k+1 , . . ., x k+m ) The popular design criteria just use the diagonal of Σ which means to voluntarily dispense of valuable information given by the kriging covariances.
The design criterion considered in this paper, is the generalized variance of the kriging covariance matrix: Definition 2. The design ξ that minimizes criterion (4) or equivalently any root of |Σ| is defined as the GV-optimal design, GV stands for Generalized (kriging) Variance.
As designs that simply minimize the kriging variance, GV-optimal designs are often space-filling but typically position more design points at the edge of the design region.GV-optimal designs depend more on the special covariance structure which is in contrast to G-and V-optimal designs in particular for small numbers of observations k (see Figure 1 for an exemplary comparison).Unfortunately the minimization of the GV-criterion is computationally demanding, since the evaluation of (4) requires the evaluation of the determinant of an (m × m)-matrix, being unfeasible for large m which is especially necessary in high dimensional design spaces as it is often the case for computer experiments.A remedy for this problem is provided by the use of incrementally assembled designs proposed here which turn out to be GV-optimal.This reduces the computational effort for the evaluation of (4) for arbitrary m to the evaluation of the determinant of a (k × k)-matrix, where k is just the design size.
The paper is organized as follows.In Section 2 we motivate our approach, exploiting the fact that the volume of the simultaneous confidence region for the prediction errors is proportional to |Σ|.Section 3 presents the basis of the main contribution of the paper, giving an update formula for the determinant of the kriging covariance matrix which should be used if we are constructing our designs incrementally.These update formulas are also necessary to speed up the computation of the design criterion (4) and may as well be used for a more efficient computation of the design Figure 1: Map of kriging variances of GV-optimal (left), G-optimal (middle) and V-optimal (right) 12-point design for the ordinary kriging setup on the unit square and Matern covariance function with κ = 2.5 and ϕ = 3.The black dots are the design points.Note that in some zones of the design region the GV-optimal design has considerable larger kriging variances than the G-and V-optimal designs.This is a natural consequence of the GV-criterion, which does not aim at minimizing the point-wise kriging variances while the G-and V-criterion do.criterion V-optimal designs.Finally, Section 4 considers the efficiency of incremental GV-optimal designs and also the efficiency of GV-optimal designs with respect to G-and V-optimal designs which is demonstrated by means of a representative simulation study in Section 5, followed by a real world example in Section 6, and Section 7 draws conclusions on the efficiency and limitations of the approach and suggests topics for future work.
2. The generalized variance and optimal designs for prediction Wilks (1932) remarked that there were ". . .statistical coefficients which have not been adequately generalized for samples from a multivariate population including the variance . . ." and introduced the generalized variance which is simply the determinant of the covariance matrix of a multivariate population.Wald (1943) used the idea of minimizing the generalized variance in his criterion for optimal designs for parameter estimation (D-optimality), i.e.
where M θ is the information matrix of the parameter estimates θ.A combination of these two ideas naturally leads to criterion (4).Shewry and Wynn (1987) introduced maximum entropy sampling (MES) where the Shannon entropy is used as a measure of information to get optimal designs for prediction.Here the design criterion is similar to (4): in the case of Gaussian response Y , MES tries to maximize the determinant of the covariance matrix of Y (ξ) which is equivalent to minimizing the determinant of the covariance matrix of (Y (x k+1 ), . . ., Y (x k+m )|Y (x 1 ), . . ., Y (x k )), i.e. the determinant of the covariance matrix of Y (•) at the unsampled locations (x k+1 , . . ., x k+m ) conditioned on Y (•) at sampled locations (x 1 , . . ., x k ).Though this criterion does not aim at minimizing the variation of the prediction errors, it is not even directly connected with a prediction method.MES it is rather a sampling method trying to absorb the maximum amount of variability into the sample, such that conditional on the sample the unsampled points have minimum variability.The method is suitable for observations on a finite closed system which is the main connection to the present work.Beyond that we connect the sampling method directly with the BLUP for the response Y as presented hereinafter.
Interestingly MES yields designs that are exactly GV-optimal in the simple kriging setup and also seem to be GV-optimal in the ordinary kriging setup.For universal kriging GV-optimal designs approach MES designs with increasing effective range, i.e. more strongly correlated response.

The best linear unbiased predictor and the corresponding kriging covariance matrix
As with MES we assume that the allowable choice of the designs ξ is from a fixed finite set of N = k + m points X = {x 1 , . . ., x k+m }.Given a k-point design ξ = (x 1 , . . ., x k ), 0 < k < N , the complementary design ξ 0 is the set X \ ξ and we get the corresponding partitioning of the response vector Y = Y T ξ Y T 0 T .We now want to simultaneously predict the response from the complementary design on the basis of the response from the design ξ.
As mentioned above we are using the model of universal kriging, i.e. our linear predictor is the linear combination of a vector of deterministic functions of the locations x ∈ ξ: ) with the first component usually being f 1 (x) = 1.The design matrix and corresponding vector of errors then are We denote the covariance matrix of ε ξ with C ξ and use analogous nomenclature for the complementary design to get The generalized least squares (GLS) estimate of the parameter vector then is β = ξ Y ξ yielding the simultaneous kriging prediction for all complementary design points as the BLUP The components of the (m × k)−weight matrix W = (ω ij ) give the kriging weights of Y (x j ) for the prediction Ŷ (x k+i ).Combining the results from above the kriging prediction errors have expectation zero and the kriging covariance matrix

Geometrical interpretation of the GV criterion
In the case of Gaussian kriging prediction errors, we get a simultaneous confidence region for the prediction errors with where χ 2 m;1−α is the (1 − α)-quantile of the χ 2 m distribution.This region is a k-dimensional ellipsoid who's volume is proportional to |Σ|, i.e.GV-optimal designs minimize the area of the prediction errors.It is evident that the usually used design criteria like G-optimality or V-optimality do not have this highly desirable property.

Comparison of MES and the GV-criterion
The criterion for a GV-optimal design or a GV-optimal increment (13) looks similar to the criterion of MES with multivariate normal distributed Y, which in our notation would be to maximize The reason for the equivalence lies in the determinant identity and the fact that the left hand side is fixed and finite (see Shewry and Wynn (1987)).
However, there are considerable differences between the GV-criterion and the MES-criterion: • The above MES-criterion just holds in the case of Gaussian Y.We do not assume multivariate normal distributed Y for our GV-criterion, i.e. the two criteria are different in the case of non-Gaussian Y.
• With GV-optimal designs we minimize Cov Ŷ0 − Y 0 in contrast to Cov (Y 0 |Y ξ ) which is minimized with MES-optimal designs.This difference just vanishes in the case of simple kriging which is shown later.
• The implementation of MES for a linear model (1) uses a Bayesian model which needs some prior distribution for the parameter vector β and particularly some prior covariance matrix Cov (β) = R −1 .The MES criterion for a linear model then turns out to be: maximize and since R does not depend on the design ξ the MES criterion for a linear model can also be formulated as: Select ξ such that Chaloner and Verdinelli (1995).We are not working in a Bayesian framework and do not assume a prior distribution or a prior covariance matrix for β, so the GV-criterion is clearly different from MES for a linear model.Despite these differences it might be interesting to compare GV-optimal designs with MES for Gaussian response just using a constant model.Here the MES-criterion is simply: choose ξ such that |C ξ | is maximized.For simple kriging this criterion is equivalent to the GV-criterion because here the kriging covariance matrix is With ordinary kriging the above equivalence is not clear: The kriging weights for ordinary kriging are where b = 1 k C −1 ξ 1 T k and 1 n is a vector of ones with length n.The kriging covariance matrix for ordinary kriging then is and the GV-criterion in this case requires the minimizing of The above criterion function is clearly different from the simple kriging setup and the following conjecture was verified by countless computations of |Σ SK | and |Σ OK | for the same randomly generated designs respectively with different covariance models, different functions f (x) and designs of different sizes.A formal proof eludes us.
Conjecture 1.Let Σ 1 SK be the GV-criterion function for design ξ 1 in the simple kriging setup and Σ 1 OK the criterion function in the ordinary kriging setup for the same design ξ 1 .Σ 2 SK and Σ 2 OK are defined analogously for another design ξ 2 .Then but the arguments of the minima are the same in both cases as can be seen in Figure 2.
In the case of universal kriging the design criterium and the optimal design is clearly different from MES even though it can be observed that with increasing effective range i.e. with higher correlated response GV-optimal designs approach MES-designs also for universal kriging.

Incrementally and decrementally constructed GV-optimal designs
In many practical situations the experiment is not stopped after a fixed number of runs but say after a certain time or when a certain budget for the runs has expired.Thus at the start of the experiment the sample size is unknown and it is not clear which optimal designs of which size to use.In these situation incremental designs should be used, starting with a first design of (small) size k and supplementing it step by step with increments of size l i ≥ 1 until experimenting has to be stopped.Here the question arises how to find optimal increments in an efficient way which will be answered in this section.The efficiency of incrementally built designs compared to GV-optimal designs will be examined in Section 4.
The use of decrementally constructed designs is not directly motivable, it arises from the attempt of detecting GV-optimal designs of size k +l in an incremental way.This may be done very efficiently with the help of Corollary 1 as will be shown below.The idea is to start with a design of minimal size k and compute an increment of size l which may be done with small computational effort.The only problem can be caused by the starting design which may contain design points which are not elements of the GV-optimal design of size k +l leading to highly efficient but still suboptimal designs.In these situations an improvement can be achieved using a decrement, i.e. omitting k 1 design points followed by another incremental step.The hope in the decremental step is to get rid of inappropriate design points which in the following incremental step are substituted by design points yielding designs closer to GV-optimality.

Update formulae for kriging weights and the kriging covariance matrix of incremental designs
For this subsection we assume that we already have a k-point design and want to add l extra design points to simultaneously predict the remaining m ≫ l non-design points.As the calculation of the design criterion is computationally demanding, we will show that the use of update formulae for the kriging weights and consequently for the kriging covariance matrix is of great computational benefit.
Furthermore these update formulae may also be used for updating Cov Ŷ (x) and Cov (Y (x)|Y (ξ)) and also for incrementally built G-and V-optimal designs.The allowable choice of the designs ξ 1 and ξ 2 with ξ 1 ∩ ξ 2 = ∅ are from a fixed finite set of N = k + l + m points X = {x 1 , . . ., x k+l+m }. ξ 1 = {x 1 , . . ., x k } is the k-point first or starting design, ξ 2 = {x k+1 , . . ., x k+l } is the l-point second design (the increment), and X \ {ξ 1 ∪ ξ 2 } is the remaining sets of non-design points with cardinality m.Thus we get the corresponding partitioning of the response vector In the first stage the design is ξ 1 and we have to predict where W 12 are the weights of Y 1 for the prediction of Y 2 and analogously W 10 the weights of Y 1 for the prediction of Y 0 .
In the second stage we add the increment and the design is now ξ 1 ∪ ξ 2 and we have to predict where W 1 are the new weights of Y 1 for the prediction of Y 0 and analogously W 2 the weights of Y 2 for the prediction of Y 0 .Emery (2009) showed that the kriging weights of the first stage (in our notation the components of W 12 and W 10 ) can be updated such that in our compact matrix notation The only problem with ( 9) is that W 2 is unknown if we add an increment to a smaller design.So, an "update" formula for the computation of the weights W 2 (in fact this is not an update because weights of Y 2 do not exist in the first stage) is essential for an efficient prediction of Y 0 in the second stage.
Emery ( 2009) also presented update formulae for the kriging variances and covariances which unfortunately are wrong in the case of l > 1, which was shown with a simple counter example by Chevalier and Ginsbourger (2012).The presented "corrected" update formulae just have not been for the kriging covariance but for Cov (Y 0 |Y ξ ).Eventually Chevalier et al. (2014) introduced correct update formulae for kriging variances and kriging covariances and also formulae for the new kriging weights of the second stage W 2 .In our notation these formulae may be summarized as follows.

(Updated) kriging weights for the second stage
Let the weights W 12 , W 10 , W 1 and W 2 be as defined in ( 7) and (8).Let further the kriging covariance matrix of the first stage in obvious notation be Then the new weights W 1 and W 2 can simply be computed with An algebraic proof for the simultaneous computation of the (m × l)-weight matrix W 2 which is fundamentally different from the one given in Chevalier et al. ( 2014) can be found in the Appendix.

Update formula for the kriging covariance matrix
Let the kriging covariance matrix of the first stage be as in (10), and the kriging covariance matrix of the second stage be Σ + 0 .Using the weights of (11) we can then update the kriging covariance matrix from the first stage to get An algebraic proof for this simultaneous update formula for the kriging covariance matrix which uses again another reasoning than Chevalier et al. (2014) can be found in the Appendix.

Efficient computation of increments for GV-optimal designs
The following theorem is of fundamental importance as it allows the efficient incremental construction of GV -optimal designs.
Theorem 2. The GV-optimal increment for the second stage given a design at the first stage then is Proof of Theorem 2: For a given design at the first stage the determinant of the corresponding kriging covariance matrix is fixed and finite and the determinant may be factored as: The GV-optimal increment is the design ξ 2 that minimizes the determinant of the kriging covariance matrix of the second stage , which is exactly the second factor of the determinant factorization above.Obviously for a fixed kriging covariance matrix of the first stage minimizing Theorem 2 is the reason for the great computational benefit of using incremental designs.It reduces the computation of the usual criterion function, which beside some matrix inversions demands the computation of the determinant of a (m × m)-matrix, to the computation of the determinant of a (l × l)-matrix with l ≪ m.This fact enables an increase of the set of points X = {x 1 , . . ., x k+l+m } to an arbitrary size N = k +l +m without raising the computational demands.
Actually, the problem of finding arg max ξ 2 |Σ 2 | is known to be NP-hard, see Ko et al. (1995), the additional demand of computing the determinant of a (m × m)-matrix would make it intractable already for moderate m > 1000.
Remark 2. Theorem 2 can even be used for an efficient computation of GV-optimal designs of a given size k + l: We start with some minimal preliminary design, i.e. a design of minimal necessary size k = p where p is the number of deterministic functions of the locations x ∈ ξ: ) used in the linear predictor.The design points can even be chosen randomly and its kriging covariance matrix is the basis for the computation of the increment ξ 2 as above.After this incremental step the design is reduced to size k 1 ≥ k.Also the k 1 design points are chosen randomly out of the incremental design of size k + l.On the basis of these k 1 points we again compute the GV-optimal increment to end up with a design of size k + l.As the computational effort is small, these decremental and incremental steps may be repeated many times.If k and l are of moderate size we may even loop systematically through all k 1 -combinations of the k + l design points.
Remark 3. By applying the above incremental step several times we may also construct highly efficient sequential designs which accounts for active learning.
We can now formulate a similar procedure for efficient computation of increments for V-optimal designs.
Corollary 3. The V-optimal increment for the second stage given a design at the first stage is Proof of Corollary 3: For a given design at the first stage the trace of the corresponding kriging covariance matrix is fixed and The V-optimal increment is the design ξ 2 that minimizes the trace of the kriging covariance matrix of the second stage tr Obviously for a fixed kriging covariance matrix of the first stage minimizing tr Corollary 3 is the reason for the computational benefit of using incremental designs.It reduces the computation of the usual criterion function which demands the matrix inversions of one (k + l) × (k + l)-and one (p × p)-matrix and the computation of the new kriging covariance matrix (10 matrix multiplications of which 5 involve matrices with m as one dimension) to the computation of the inverse of a (l × l)-matrix and 2 matrix multiplications of which only 1 involves a (l × m)-matrix.
As l ≪ m this reduces the computational effort to roughly one third.Just as Theorem 2, Corollary 3 can be used for an efficient computation of designs close to Voptimality of a given size k + l.The computational benefit of Corollary 3 cannot be compared to the improvement of Theorem 2, as we have to limit the number of incremental and decremental steps here and we thus have no guarantee to end up with the V-optimal design.

GV-optimal designs in design spaces dense in R d
In subsection (2.1) we state that the allowable choice of the designs ξ is from a fixed finite set of N = k + m points X = {x 1 , . . ., x k+m }.Usually k ≪ m and going into details of Corollary 1, Remark 2 we may even extend the number of non-design-points m to any (integer) size.
The reason for this remarkable possibility is that for the computation of the GV-optimal increment ξ 2 of a given design ξ 1 we do not even need the kriging covariance matrix of the first stage which would be of enormous dimension (l + m) × (l + m), we just need the (l × l)-block Σ 2 which is The dimensions of all of the above matrices are just k or l.W 12 are the weights of Y 1 for the prediction of Y 2 in the first stage: The dimensions of all these matrices are only k, l and p, i.e. for the computation of Σ 2 we just need small matrices independent of the number m.
The GV-criterion function is the determinant of the kriging covariance matrix, if in the first stage this determinant is D, then the criterion function of the incremental design is D |Σ 2 | .In the decremental step we remove l 1 = k + l − k 1 design points from the incremental design, the according (l 1 × l 1 )-block of the kriging covariance matrix of the first stage is Σ * 2 and may be computed as shown above.
The determinant of the kriging covariance matrix after the decremental step then is D |Σ 2 | •|Σ * 2 |.So, not knowing the value of the GV-criterion during the search for the GV-optimal design is not crucial as it suffices to know the Σ 2 -blocks according to the increments and decrements to minimize the GV-criterion.Thus in principle we can make the grid as dense as desired, as long as the number of points is finite.There is reasonable hope that the described method can be generalized to continuous design spaces, which we defer to future research.

Efficiency of GV-optimal designs
As mentioned above it turns out that GV-optimal designs are highly efficient with other design criteria which will be discussed in subsection (4.1).
Additionally Corollary 1 allows a simple, fast and computationally very efficient calculation of incrementally constructed designs that are close to be GV-optimal.The efficiency of these incremental designs will be discussed in subsection (5.1).

Efficiency with respect to other design criteria
Traditional criteria for optimal designs for prediction are usually concerned with the variance of predictions, i.e. we could also title this subsection with "Efficiency with respect to variance-based criteria".Here we compare our GV-optimal designs with G-and V-optimal designs with the help of the relative efficiency, a very common concept in comparing designs, see eg.López-Fidalgo (2023), p.17.
Let Φ G , Φ V and Φ GV be the criterion function for G-, V-and GV-optimality respectively and ξ G , ξ V and ξ GV be the G-, V-and GV-optimal designs of the same size for the prediction of the same number of points.The GV-optimal design ξ GV minimizes Φ GV , as the other optimal designs minimize their corresponding design criteria.Then e.g. the relative G-efficiency of the V-optimal design is We always have E • (•, •) <= 1 and the relative efficiency of a design ξ gives the factor the criterion function of ξ may be decreased if we switch to the optimal design.These relative efficiencies are scale invariant though the effect of scaling is different for the GV-criterion function on the one and the Gand V-criterion functions on the other hand.The kriging covariance matrix of scaled responses s • Y is s 2 Σ which affects G-and V-criterion functions the same: max diag (s 2 Σ) = s 2 max diag (Σ) and tr (s 2 Σ) = s 2 tr (Σ).The GV-criterion can be made insensitive to scaling by applying k |Σ| instead without changing the designs.The relative efficiencies are not affected by this scaling.I.e., if Σ GV , Σ G and Σ V are the kriging covariance matrices for the originally unscaled data of the GV-, G-and V-optimal designs, then the GV-efficiency of the G-and V-optimal designs respectively for scaled data are i.e. arbitrary scaling does not change the relative GV-efficiencies.The same is true also for relative G-and V-efficiency.

A representative example
Let us demonstrate the typical relative efficiencies on the basis of the following settings, computations of many other differently adjusted models and designs yield similar results.
As can be seen in tables 1 and 2, GV-optimal designs are reasonably efficient with respect to the G-and V-criterion.Here for every covariance parameter combination the GV-, G-and V-optimal designs were determined and then for each optimal design the relative efficiencies with respect to the other design criteria were computed.This was here done for a linear trend function and for designs of size 6, 9 and 12 respectively.Finally, the relative efficiencies of the optimal designs were averaged over all 54 covariance parameter combinations.The lines of the tables correspond to GV-, G-and V-optimal designs, and e.g. the mean relative efficiency of 9-point GV-optimal designs with respect to the G-criterion is 0.9707 which means that the maximum kriging variance of G-optimal designs is on average only 97% of the maximum kriging variance of GV-optimal designs.linear 6 points 9 points 12 points trend

Efficiency of incrementally assembled designs
As already mentioned incrementally assembled designs are very efficient with respect to the GVcriterion.In the following discussion we will always start with some k-point design ξ k adding a single increment of size l.The result will only be the GV-optimal design of size (k + l) if we start in the first stage with a k-point design ξ k ⊂ ξ k+l which is very improbable if we do not utilize additional knowledge.Usually we also will not end up in the GV-optimal (k + l)-point design if we start with the GV-optimal k-point design, but the result will be very close to GV-optimality.How close the incrementally constructed design is to GV-optimality depends on the choice of the k-point starting design.Of course we may take advantage of prior knowledge about properties of GV-optimal designs, e.g. if the design region is the unit square as in our simulation examples, we know that the GVoptimal design for a linear or quadratic trend will always have design points in the corners and the edges of the design region.Choosing such plausible points for the k-point starting design will almost always yield GV-optimal (k + l)-point designs that are constructed with a single incremental step.
Here we followed two ideas: • start with a GV-optimal design of (small) size k; • start with a plausible design of (small) size k, i.e. with design points that most likely are elements of GV-optimal designs of arbitrary size.
It turns out that both ideas yield very efficient designs especially if the increment (the number of additional design points) is not too small.
To exemplify this efficiency we again used Matern covariance models with 54 different parameter combinations and a quadratic trend.In the first simulation series we started with a plausible 6-point design, i.e. 4 design points at the corners and 2 points on opposite margins of the unit square.Then we added the optimal increment of 6 design points as described above.The mean GV efficiency of these incremental designs was 0.9911, the median efficiency was even 100%.
In the second series of simulations we started with the 7-point GV-optimal designs for each parameter combination respectively and added the optimal 5-point increment.Here the mean GV efficiency was 0.9862 and the median efficiency again 100% indicating a satisfactory performance of this simple incremental method.
We applied the same concept to linear trends as well.The 4 corners of the design region were chosen as plausible starting design, and the increment of size 8 always yielded the GV-optimal 12point design.Starting with 6-point GV-optimal designs for each parameter combination respectively and adding the optimal 6-point increment resulted in a mean GV-efficiency of 0.9881, again the median efficiency was 100%.
Note that the above mean and median efficiencies are just for a design built with a single incremental step.Of course we may always append a few decremental-incremental steps to guaranty GV-optimality.This approach is analyzed in the next section.

Efficiency of incrementally-decrementally assembled designs
The above described method of systematically discarding design-points after each incremental step has an efficiency of 100%.For each of the 54 combinations of the κ and ϕ parameters of the Matern covariance model (see above) we started 1000 times with a random design and ended with the GV-optimal design every time just with 26 exceptions (where a design with efficiency 99.9% was found instead of the optimal design).It turned out that theses exceptions were all for designs with three special parameter combinations of κ and ϕ which seem to be adverse for finding the GV-optimal design on the chosen grid.For these three parameter combinations also the average number of required computations of the criterion function was incomparably higher then for other parameter values (see Table 3).The reason was obviously that the design points were limited to the unfavourable (17 × 17) grid.Changing to a (33 × 33) grid solved this problem.With the finer grid we started 200 times for each of the 54 parameter combinations of κ and ϕ with a random design and found the optimal designs without exception.Also the number of computations of the criterion function till convergence was distributed more uniform than with the (17 × 17) grid.Though there were some parameter combinations which needed clearly more function calls, for all these cases the second best design always was very efficient and the search algorithm now and then got stuck at these designs before finding the optimum.Table 3: Median number of computations of the criterion function until convergence to the GV-optimal design.

Speed of convergence to GV-, G-and V-optimal designs
The speed of convergence was measured in absolute time and in the number of required computations of the criterion function.
In Table 3 we can see the median number of computations of the criterion function (in this case the determinant of the (6 × 6)-matrix Σ 2 ) needed to find the GV-optimal design.The overall median number of calls of the criterion function was 17222, the computation time for finding 54.000 times the GV-optimal designs was 50.14 hours; i.e. 3.34 seconds per optimal design.
The distributions of the number of criterion function evaluations turned out to be positively skewed, the cdf of such distributions for selected parameter combinations of κ and ϕ is depicted in Figure 3.There were 3 parameter combinations where the corresponding distributions of function calls were striking.Discarding these extreme distributions reduced the average computation time for one GV-optimal design to 2.18 sec.Changing to a grid twice as fine as the original (17 × 17) points to which the design was limited solved the problem (see above).The overall median number of calls of the criterion function increased to 24877 with the (33 × 33) grid which was caused by the halved step size in the neighbourhood search at the finer grid.Here continuous optimization algorithms with variable step size promise an improvement (see Section 3.3).
The GV-optimal designs are found incomparably faster than corresponding G-and V-optimal designs even if the allowable choice of designs is from a moderate number of points.
Searching for G-optimal designs could only be tried 10 times for each of the 54 combinations of the κ and ϕ parameters and 368 of these 540 tries failed, because we had to stop the search algorithm   (a combination of neighboring point exchanges and simulated annealing, same algorithm was used to find the GV-optimal designs) after 500 iterations because of the huge expenditure of time.The reason for that was partly the much larger computational effort but also a much slower convergence to the optimum, i.e. the criterion function (max diag Σ 0 − Σ T 20 Σ −1 2 Σ 20 ) had to be called much more frequently than in the search for GV-optimal designs.The mean G-efficiencies of the 540 found designs was 0.991, for κ and ϕ parameter combinations corresponding to high correlated data the mean G-efficiencies of the found designs were considerably smaller (with a minimum of 0.927 for κ = 2.5 and ϕ = 5).
In Table 4 the median number of calls of the criterion function for each of the 54 combinations of the κ and ϕ parameters is reported.The overall median number of calls of the criterion function was 5.101.156(296.2 times as often as for GV-optimal designs), the computation time for finding 540 times the G-optimal designs was 761.5 hours, i.e. 1.41 hours per optimal design which is 1520 times as long as for one GV-optimal design.
V-optimal designs are somehow found easier than G-optimal designs (because we may apply Corollary 2).We managed to search the V-optimal design 250 times for each of the 54 combinations of the κ and ϕ parameters.In table 5 we can see the median number of computations of the criterion function (in this case tr (Σ 2 ) + tr Σ −1 2 Σ 20 Σ T 20 ) needed to find the V-optimal design.Also here the average number of calls of the criterion function was clearly larger than for the GV-optimal design.The overall median number of calls of the criterion function was 108928.5 (6.3 times as often as for GV-optimal designs), the computation time for finding 13.500 times the V-optimal designs was 240.7 hours, i.e. 64.2 seconds per optimal design which is 20 times as long as for one GV-optimal design.
Of the 13.500 tries to find the V-optimal design 577 failed, that is 4.3% (almost exactly 100 times more than for GV-optimal designs).This has also an impact on the cdf of the positively skewed distribution of the number of calls of the criterion functions until the optimal designs were found (Figure 4).As with G-optimal designs we stopped the search algorithm after 500 iterations, the cdf's for parameter combinations where the optimal designs were not found within this maximal number of iterations are depicted as colored lines.
Figure 4: Empirical cdf of the number of computations of the criterion function till convergence to the V-optimal design (black line) and for parameter combinations where the optimal design was not always found within a maximal number of tries (colored lines) -mind the log-scale.

Real illustrative example: temperature prediction in Upper Austrian municipalities
The province of Upper Austria is partitioned in 438 rural and urban municipalities with considerable topographical differences ranging from lowlands in the center and hill country in the north to high-altitude mountains in the south.As temperature and its spatial variation is strongly influenced by the topography the simultaneous prediction of temperatures at all principle locations of the 438 municipalities is challenging.
Currently there exist 36 meteorological stations in Upper Austria that may be taken as data source for temperature prediction in the 438 municipalities.A natural question is whether the current network can be improved by relocation of the station and/or we can even reduce the size of the network without loss of accuracy.
We model the expected monthly mean temperatures t ij with the elevation of the measurement location el i as external drift which is in line with Hudson and Wackernagel (1994): where j indicates the month and i the location of the measurements.We further use an anisotropic Matérn covariance model to describe the spatial interdependencies of temperatures measured in the same time period.The covariance parameters have been estimated with the likfit function of the R package GeoR (Ribeiro Jr. and Diggle (2001)).
The learning data for parameter estimation were the daily mean temperatures of all meteorological stations in Upper Austria in the period from 2000-01-01 until 2023-10-25.The data are publicly available at the GeoSphere Austria Data Hub (2023).The coordinates and elevations of the 438 municipalities are also publicly available at the DORIS webOffice (https://www.doris.at/) The parameter estimates confirm the environmental temperature lapse rate of ∼ 6.5°C/km (International Civil Aviation Organization (1993), Thompson (1998)) and are similar for all months except the winter period when the phenomenon of temperature inversion (National Oceanic and Atmospheric Administration's (2023)) may be observed frequently.
Here the showcase results for the month June are presented, other months are comparable.The design region is the set of 438 locations corresponding to the principle localities of the 438 Upper Austrian administrative municipalities, actually 36 of these locations are the base of meteorological stations (https://bitly.ws/ZFpC).We want to evaluate the prediction quality of this actual meteorological network by means of the GV-criterion and compare it with a virtual network positioned at the locations of the GV-optimal design of the same or a reduced size.I.e., initially we have k = 36 and m = 438 − 36 = 402 locations for simultaneous prediction.
We further reduced the size step by step down to k = 32 until which we yielded the same accuracy according to the GV-criterion as from the the initial network.This resulting GV-optimal design as well as the actual meteorological network and all 438 locations are displayed in Figure 5. Remarkably 7 locations of the optimal design are already base of a meteorological station and other 7 locations of the optimal design are within a range of 5km from an actual station although to our knowledge no statistical analysis was involved in determining the current network.While the reduction of just 4 stations from 36 down to 32 may seem as a disappointment one must not forget that the number of prediction locations had to increase from 402 to 406 making the task more difficult.

Conclusions
With GV -optimality we have introduced a novel design criterion for simultaneous kriging prediction, which considers the whole prediction covariance matrix.As was shown by the real-world example there are indeed practical problems requiring simultaneous rather than individual prediction.In such situations, the presented new criterion is a natural answer and more adequate and useful.
In terms of robustness, it has been demonstrated that GV-optimal designs exhibit considerable efficiency compared to designs optimized based on other criteria.The criterion function is notably smoother compared to other criteria, meaning that slight modifications to the design do not lead to significant alterations in the criterion function.Interestingly, this might also be the main reason that GV-optimal designs are found much faster than designs optimal with respect to other criteria (which is subject to actual and future research).
Furthermore we have shown that efficient incremental construction methods are available, which makes the criterion particularly attractive for big data and higher dimensional contexts.For instance it lends itself naturally combinable with local kriging techniques such as Gramacy and Haaland (2016).These and other extensions will be subject of future research.

Figure 2 :
Figure 2: GV-criterion function for simple kriging |Σ SK | against GV-criterion function for ordinary kriging |Σ OK | for different designs.The minima of both criteria seem to be the same.

Figure 5 :
Figure5: Map of the 438 municipalities (dots), the current 36 meteorological stations (crosses), the locations of the 32-point GV-optimal design (red circles) for a meteorological network for simultaneous prediction of temperatures in the 438 municipalities of Upper Austria.

Table 1 :
Average relative efficiencies for linear trend over 54 different parameter combinations of Matern covariance models.

Table 2 :
Average relative efficiencies for quadratic trend over 54 different parameter combinations of Matern covariance models.

Table 4 :
Median number of computations of the criterion function until convergence to the G-optimal design.