A lot of computational models applied in the business and industry, using multidimensional data, obtained by measuring, sometimes require to fulfil expensive simulations. In these cases, the kriging method is proposed to increase the efficiency of solving modelling, prediction and optimization problems (Jones, 2001; Forrester and Keane, 2009; Xiao et al., 2018). Since the kriging method was developed as a consequence of its accurate predictions (Krige, 1951), it became one of the most promising numerical approaches in the various fields of engineering design, spatial statistics, experimental design, and so on (Kwon and Choi, 2015; Bhosekar and Ierapetritou, 2018; Carpio et al., 2017). In general, kriging is a method of interpolation for which the interpolated values are modelled by a multivariate Gaussian process, governed by prior covariances. The basic idea of kriging is to predict the value of a response function at a given point by computing a weighted average of the known values of this function in the neighbourhood of the point. Thus, kriging gives a way of anticipating, with some probability, a result associated with values of the parameters that have never been met before using the existing information (e.g. the experimental measurements).

The first step in kriging is to create a stochastic model that likely best describes the set of observed data. The assumption is usually made that the sample values are sufficiently homogeneous, as well as the correlation between two random points solely depends on the distance between them, but it is independent of their location. In this way, the stochastic homogeneous and isotropic model for kriging is constructed using the correlation functions that have to obey the conditions of positive definiteness according to Bochner’s theorem (Abrahamsen, 1997). Following this approach, many models in spatial statistical research and applications are developed using the Whittle-Matern correlation family (Abrahamsen, 1997; Guttorp and Gneiting, 2005), Gneiting correlation class (Guttorp and Gneiting, 2005), etc. In this paper, the application of Euclidean distance matrices (EDM) to multivariate data modelling and kriging is considered, doing it without the standard way based on positive definite correlation functions.

EDM’s have received increased attention due to their applications in recently active fields of research, such as molecular conformation in bioinformatics, dimensionality reduction in machine learning and statistics, semidefinite programming, wireless sensor network localization (Pham et al., 2008; Ghiasi et al., 2018; Koo et al., 2018; Qian et al., 2008). However, in recent years square EDM’s are mostly studied in literature (Schoenberg, 1935; Gower, 1984; Weinberger et al., 2004), although the Euclidean geometry in the real world is dealing with a square root, i.e. the fractional degree of distance squares. Moreover, fractional degrees of distance squares appear in some applications as well, however many fundamental issues about their properties remain unstudied. Therefore, this work focuses on the geometry of Fractional Euclidean Distance Matrices (FEDM) and their application in the multidimensional data analysis using kriging in the context of computer modelling and solving of practical examples. Sakalauskas (2013) considered a random field with a positive definite correlation function tending its scale parameter to zero as well as infinitely increasing the variance of field, and in this way has derived that the limiting kriging method relates with inverse FEDM’s (see, also Remark 1 below). Usually EDM’s are studied using the Gram matrix of the observed data set, called a kernel matrix (Gower, 1984, etc.). Using the kernel matrix method (Pozniak and Sakalauskas, 2017) has shown FEDM with a strongly fractional degree of distance squares to be nonsingular, and, consequently, the derived kriging method in can be used correctly. In this paper, the multinormal data model and the kriging method are developed taking covariances proportional to the elements of kernel matrix, and, thus, avoiding the standard way, based on positive definite correlation functions. Once constructed, the kriging model can be used for data modelling, prediction, and optimization. The objective is to discuss and demonstrate the applicability of the kriging technique, based on FEDM, and to verify the computational efficiency of the developed approach.

Let us assume that a data set X = ( x 1 , x 2 , … , x K ) of K vectors x i ∈ R d , 1 ⩽ i ⩽ K, d ⩾ 1 be given, at which the values of some response function Y = ( y 1 , y 2 , … , y K ) T are used, obtained by physical measurement, computer simulation, etc. Denote the K × K matrix of fractional degrees δ of Euclidean distance squares among pairs of vectors by A = [ ( ‖ x i − x j ‖ ) δ ] i , j = 1 K , where ‖ x i − x j ‖ = ( x i − x j ) T · ( x i − x j ), x i ∈ R d , 1 ⩽ i , j ⩽ K, 0 ⩽ δ ⩽ 1. Note that δ = 1 / 2 in a special case of usual Euclidean distances. Recently the properties of matrices of Euclidean distance squares are studied most, i.e. as δ = 1 (Schoenberg, 1935; Gower, 1984; Weinberger et al., 2004). Usually the main tool for EDM study is the Gram matrix of data set X, called a kernel matrix. Thus, following this approach, the matrix (1) F = − ( I − E · s T ) · A · ( I − s · E T ) is introduced, called hereinafter a kernel matrix, where I denotes K × K unit matrix, E denotes K-dimensional vector-column of units, s ∈ R K is a certain vector-column.

The proof can be seen in Pozniak and Sakalauskas (2017), Theorem 1. Note that, if δ = 1, then the rank of kernel matrix can be less K − 1.

The proposition implied by an easily verified equality of idempotence as ( I − E · s 1 T ) · ( I − E · s 2 T ) = ( I − E · s 1 T ) .

Let us consider FEDM geometrical properties more.

Hence, Corollary 1 shows us the geometrical sense of vector s, that displays the origin of coordinates to the point z c in the embedded space. Therefore, if s T · E = 1, it is reasonable to call s a centering vector.

The next corollary follows from (1).

Let us consider the inverse and determinant of FEDM more in detail. Denote by ( x 0 , X ) the set obtained adding some point x 0 ∈ R d to the observation set X, as well as the respective vector of fractional degrees of Euclidean distance squares between x 0 and X by a = ( ‖ x 1 − x 0 ‖ δ , ‖ x 2 − x 0 ‖ δ , … , ‖ x K − x 0 ‖ δ ) T . Now, define the extended kernel matrix, respective to FEDM A ˜ = 0 a T a A of the set ( x 0 , X ), calculated at the centering vector s ˜ T = ( s ′ , s T ) and presented in the block matrix: (2) F ˜ = ν f T f F ≡ − 1 − s ′ − s T − s ′ · E S · 0 a T a A · 1 − s ′ − s ′ · E T − s S T , where, due to (1): (3) S = I − E · s T , (4) F = s ′ · E · a T · S T + s ′ · S · a · E T − S · A · S T , (5) f = S · A · s − s ′ · E · a T · s − S · a · ( 1 − s ′ ) , (6) ν = 2 · ( 1 − s ′ ) · s T · a − s T · A · s .

Then next theorem relates the inverse and determinant of FEDM with that of the kernel matrix.

The proof of Theorem 1 and Corollary 5 is given in Appendix. It is of interest that the right part of conclusion of Corollary 5 is invariant to the choice of the centering vector s and point x 0 .

Statistical approach for modelling the response function by a random field has the origin in the works of Kushner (1964), Žilinskas (1985), and Mockus (1989). In this section, a probabilistic model of response surface (scalar function) is developed, the values of which are collected by observation, physical measurements or computer simulation, etc. Due to the fact that there are no further data except the measurement performance at the experimental points, the surface which represents the objective response function can be designed as a homogeneous Gaussian random field (GRF) Z ( x , ω ), which for each point in the variable space x ∈ R d , is a measurable function of random event ω ∈ ( Ω , Σ , P ) in some probability space (Mockus, 1989; Jones, 2001; Sakalauskas, 2013). To model the observed data in a probabilistic way, one has to define the probability distribution of response surface values Z ( x , ω ) = ( Z ( x 1 , ω ) , Z ( x 2 , ω ) , … , Z ( x K , ω ) ) in the given set of points X = ( x 1 , x 2 , … , x K ). Consider this probability distribution with the constant mean vector: (12) E Z ( x , ω ) = μ · E , and covariance matrix (13) E ( Z ( x , ω ) − μ · E ) · ( Z ( x , ω ) − μ · E ) T = β 2 · F , where μ and β are parameters, β > 00$]]>, and F is a positively defined matrix.

Taking into account geometrical properties of FEDM, explored in the previous section, choose matrix F as a lower-right submatrix of kernel matrix (2), calculated using some point x 0 and centering vector ( s ′ , s ) such that s ′ + E T · s = 1, s ′ ≠ 0. Note that matrix F is positively defined according to Proposition 1.

Thus, the multivariate Gaussian density is written down as well: (14) p x 0 , X ( Y ) = exp [ − ( Y − μ · E ) T · F − 1 · ( Y − μ · E ) 2 · β 2 ] ( 2 π ) d / 2 · β d · | F | 1 / 2 .

However, as one can easily see that this density, using denotations of the previous section and (7)–(11), can be rewritten in an equivalent manner: (15) p x 0 , X ( Y ) = exp [ − 1 2 · β 2 · Y T ( A − 1 − A − 1 · E · E T · A − 1 E T · A − 1 · E ) · Y − ( Y T · q − μ ) 2 2 · β 2 · D ] ( 2 π ) K / 2 · β K · ( ( − 1 ) K + 1 · | A | · D · E T · A − 1 · E ) 1 / 2 .

In order to describe the random field by its finite-dimensional (cumulative) distributions, latter ones have to obey the consistency conditions of symmetry and compatibility (Abrahamsen, 1997). This is done by the next theorem under the appropriate choice of centering vector. Theorem 2.

The random Gaussian field Z ( x , ω ) exists in some probability space ( Ω , Σ , P ), having the functions (16) p x 0 , X K ( Y K ) = exp [ − 1 2 · β 2 · Y K T ( A K − 1 − A K − 1 · E K · E K T · A K − 1 B K ) · Y K − ( Y K T · A K − 1 · ( a K + E K · C K B K ) − μ ) 2 2 · β 2 · ( D K − C K 2 B K ) ] ( 2 π ) K / 2 · β K · ( ( − 1 ) K + 1 · | A K | · ( D K · a K · B K − C K 2 ) ) 1 / 2 as its finite-dimensional densities of distribution of Y K = ( y 1 , y 2 , … , y K ) T , y i = Z ( x i , ω ), where B K = E K T · A K − 1 · E K , C K = 1 − E K T · A K − 1 · a K , D K = a K T · A K − 1 · a K , X K = ( x 1 , x 2 , … , x K ) is a sequence of mutually disjoint points x i ∈ R d , and disjoint with x 0 ∈ R d as well, 1 ⩽ i ⩽ K, A K = [ ( ‖ x i − x j ‖ ) δ ] i , j = 1 K , a K = ( ‖ x 1 − x 0 ‖ δ , ‖ x 2 − x 0 ‖ δ , … , ‖ x K − x 0 ‖ δ ) T , E K = ( 1 , 1 , … , 1 ) T , E K ∈ R K , K = 2 , 3 , … , μ and β are parameters, β > 00$]]>.

The stochastic model of kriging should incorporate uncertainty about quantities in unobserved points and to quantify the uncertainty associated with the kriging estimator. Kriging gives us a way of anticipating, with some probability, a result associated with values of the parameters that have never been met before, or have been lost, to “store” the existing information (the experimental measurements), and propagate it to any situation where no measurement has been made. As it is unknown which of all function variables will be preponderant, consider them as equivalent, and thus calculate a distance between the measurement points, which now are symmetric with respect to the miscellaneous variables. Suppose that we have to predict the value of response surface y at some point x ∈ R d , if the set X = ( x 1 , x 2 , … , x K ) of observed mutually disjoint vectors x i ∈ R d , 1 ⩽ i ⩽ K, K > 11$]]>, d ⩾ 1, is fixed and the data of measurement Y = ( y 1 , y 2 , … , y K ) T of the response surface at points of X be known. Let the data be modelled by the probabilistic Gaussian model with constant mean (12) and the covariance matrix of type (13). Then, decompose the covariance matrix β ˆ 2 · F ‾ of vector Y ‾ = ( Y , y ) in to the covariance matrix β ˆ 2 · F of the vector Y, to vector β ˆ 2 · f of covariances between vector Y and y, and to variance of y, denoted by ν: (17) β 2 · F ‾ = β 2 · F f f ν .

According to the introduction of kriging (see Jones, 2001) it is defined by the kriging predictor (18) y ( x ) = Y T · F − 1 · ( f + E · ( 1 − E T · F − 1 · f ) E T · F − 1 · E ) , a posteriori variance of the predictor (19) β 2 ( x ) = β ˆ 2 · ( ν − f T · F − 1 · f + ( 1 − E T · F − 1 · f ) 2 E T · F − 1 · E ) , and the MLE parameter of variance (20) β ˆ 2 = 1 K · ( Y T · F − 1 · Y − ( Y T · F − 1 · E ) 2 E T · F − 1 · E ) .

Let us consider the kriging method based on the probabilistic data model, derived in the previous section.

Proof is given in Appendix.

Thus, kriging predictor (21) turns out to a linear extrapolator y ( x ) = Y T · u ( x ), with the extrapolation weights: (24) u ( x ) = A − 1 · ( a + E · ( 1 − E T · A − 1 · a ) E T · A − 1 · E ) , and variance (22) helps us the measure of accuracy of extrapolation.

The next property of the kriging extrapolator easily follows from (24).

Besides, it is easy to make sure that kriging predictor (13) coincides with the values of response at measured points. Remark 1.

The kriging model has also been derived by Sakalauskas (2013) studying the GRF with the positive defined covariance function of the shape γ δ · β 2 α · ( 1 − α · ( r γ ) δ + O ( r γ ) 2 · δ ) and tending the parameter γ to infinity.

The well-known approach for data interpolation is presented by the Shepard method (Shepard, 1968): y Shepard ( x ) = Y T · u ( x ) , where (26) u ( x ) = w ( x ) E T · w ( x ) , if | x i − x | ≠ 0 ∀ i , y i , if | x i − x | = 0 , the weights are chosen in the following way: (27) w ( x ) = ( 1 | x 1 − x | δ , 1 | x 2 − x | δ , … , 1 | x K − x | δ ) .

It is easy to see that the Shepard extrapolator also satisfies the condition E T · u ( x ) = 1.

A set of analytic test functions was chosen aiming to compare the developed kriging predictor with the Shepard extrapolator (Table 1).

At the beginning of the experiments, the behaviour of response of the system is generally unknown, so different analytic functions should be designed. In total six types of functions are considered (Kwon and Choi, 2015). The first test function is a Branin function (fourth-order polynomials) and it shows a dominant second-order trend. This function has an extremely complex and highly nonlinear behaviour. The second test function is a Linear function composed of polynomials and trigonometric functions, and it shows a strong first-order trend. The Rosenbrook function shows both the first and second-order trends.

Each test function is shown in Table 1. The nonlinear Haupt function is composed of trigonometric functions. The test function results indicate a second-order trend of both independent variables x and y.

Each test function is plotted (see Figs. 1–6).

Model parameters β 2 and δ have been estimated by the maximum likelihood method (23) and (25). The computer simulation experiment has been performed generating N = 200 samples of K = 20 randomly simulated points for each test function in its domain.

In the prediction problem, the value of the Euclidean fraction degree δ could be chosen so as to get the best model results. An elegant alternative is to calculate the optimal δ value according to (25).

Mean and standard deviation of variance parameter β 2 and Euclidean fraction degree δ are given in Table 2, histograms are presented in Figs. 7–18 ( K = 20, N = 200).

To compare the accuracy of the model the True Error Criterion is introduced (Kwon and Choi, 2015). It is defined as follows: (34) TE ( y ( x ) ) = ∑ ∑ ( y ( x ) − z ( x , y ) ) 2 , (35) TE ( y Shepard ( x ) ) = ∑ ∑ ( y S h e p a r d ( x ) − z ( x , y ) ) 2 .

The results are summarized in Table 3, and the best estimation of the true error are marked in bold for each test function.

True error comparisons results are as follows K = 20, N = 200.

Therefore, the kriging predictor, being a posterior expected value of Gaussian random field, presents itself as an efficient extrapolator of scattered data, and, in turn, the variance of kriging predictor is an efficient measure of prediction or extrapolation error.

Inadequate treatment of surface wastewater may considerably impair the quality of water. With increasing urbanization, intensification of car traffic and increasing area of impervious surfaces, the pollution of surface water and a negative impact on the aquatic environment are rising as well. One of surface wastewater treatment technologies, capable of reducing suspended solids, heavy metals and other pollutants, is surface wastewater treatment filters (Meškauskaitė and Marčiulaitienė, 2016; Meškauskaitė, 2017). The effectiveness of filters, filled with construction waste and biocarbon, was analysed using the kriging method with distance matrices. The developed method allows modelling filter characteristics with different filler ratios, based on the previous experimental studies of filters.

Let x = ( x 1 , x 2 , x 3 , x 4 ) T be a vector of filler ratios. Denote different filter characteristics: Y i j , 1 ⩽ i ⩽ m, 1 ⩽ j ⩽ K, where K = 4 is the number of experiments, m – is the number of filter characteristics, that describe capability to treat different wastes, depending on filler ratios of the filter.

Filter fillers:

1 – Quartz sand;

2 – Shredded autoclaved aerated concrete (66.7%) and stone wool (33.3%);

3 – Shredded autoclaved aerated concrete (33.3%) and biocarbon (66.7%);

4 – Shredded autoclaved aerated concrete (33.3%), biocarbon (33.3%) and stone wool (33.3%).

These filters with different fillers are designed for treatment of the main pollutants of surface wastewater: zinc (Zn), copper (Cu). Then the experiment matrix of filler proportions was chosen as follows: X = 1 0 0 0 0 0.667 0 0.333 0 0 0.333 0.667 0 0.333 0.333 0.333 .

After the experiment, the filtration characteristics have been presented by the measurement matrix Y = ( Y 1 , Y 2 ) T (see Table 4).

Filters data extrapolation, based on experiments by the kriging method with distance matrices (20), (21), has been perfomed, where FEDM A = [ ( ( x i , 1 − x j , 1 ) 2 + ( x i , 2 − x j , 2 ) 2 + ( x i , 3 − x j , 3 ) 2 + ( x i , 4 − x j , 4 ) 2 ) 1 / 2 ] i , j = 1 K and the vector of distances from the selected extrapolation point x to the experimental points for i = 1 , 2 , 3 , 4 is τ i ( x ) = ( ( X i , 1 − x 1 ) 2 + ( X i , 2 − x 2 ) 2 + ( X i , 3 − x 3 ) 2 + ( X i , 4 − x 4 ) 2 ) 1 / 2 . Visualization of the cleaning efficiency (in %) for the fillers proportion (BC = 1–QS–SHAAC–SW) using the kriging approach proposed, is given in Figs. 19–22.