Kriging Model

The Kriging model has its roots in the field of geostatistics—a hybrid discipline of mining, engineering, geology, mathematics, and statistics (Cressie, 1993)—and is useful in predicting temporally and spatially correlated data. The Kriging model is named after D. G. Krige, a South African mining engineer who, in the 1950s, developed empirical methods for determining true ore grade distributions from distributions based on sample ore grades (Matheron, 1963). Several texts exist that describe the Kriging model and its usefulness for predicting spatially correlated data (Cressie, 1993) and mining (Journel and Huijbregts, 1978). Kriging meta models are extremely flexible because you can choose from a wide range of correlation functions to build the meta model.

In addition, depending on the correlation function that you choose, the meta model can either “honor the data,” providing an exact interpolation of the data, or “smooth the data,” providing an inexact interpolation (Cressie, 1993).

The following commonly used Kriging correlation functions are available in Isight:

Exponential
Gaussian
Matern Linear
Matern Cubic

Kriging postulates a combination of a polynomial model and departures of the following form:

y (x) = f (x) + Z (x)

where $y (x)$ is the unknown function of interest, $f (x)$ is a known polynomial function of $x$ called the trend, and $Z (x)$ is the realization of a stochastic process with mean zero, variance $σ^{2}$ , and nonzero covariance. The $f (x)$ term in the previous equation is similar to the polynomial model in a response surface, providing a “global” model of the design space. In many cases $f (x)$ is taken to be a constant term $β_{0}$ ; the Isight implementation assumes a constant $f (x)$ term.

While $f (x)$ “globally” approximates the design space, $Z (x)$ creates “localized” deviations so the Kriging model interpolates the $n_{s}$ sampled data points. The covariance matrix of $Z (x)$ that dictates the local deviations is

C o v [Z (x_{i}), Z (x_{j})] = σ^{2} R ([R (x_{i}, x_{j})])

where $R$ is the correlation matrix, and $R (x_{i}, x_{j})$ is the correlation function between any two of the $n_{s}$ sampled data points $x_{i}$ and $x_{j}$ . $R$ is a $[n_{s} \times n_{s}]$ symmetric, positive definite matrix with ones along the diagonal.

Many different correlation functions exist. The following correlation functions are provided within Isight:

Name	Correlation
Exponential	$c o r r (X_{i}, X_{j}) = \prod e^{{- θ}_{k} \| X_{i k} - X_{j k} \|} .$
Gaussian	$c o r r (X_{i}, X_{j}) = \prod e^{{- θ}_{k} {\| X_{i k} - X_{j k} \|}^{2}}$
Matern Linear	$c o r r (X_{i}, X_{j}) = \prod (1 + θ_{k} \| X_{i k} - X_{j k} \|) e^{{- θ}_{k} \| X_{i k} - X_{j k} \|}$
Matern Cubic	$c o r r (X_{i}, X_{j}) = \prod (\frac{1 + θ_{k} \| X_{i k} - X_{j k} \| + \frac{1}{2} θ_{k}^{2} \| X_{i k} - X_{j k} \|}{3}) e^{{- θ}_{k} \| X_{i k} - X_{j k} \|}$

In the above table $n_{d v}$ is the number of design variables and $θ_{k}$ are the unknown correlation parameters used to fit the model.

Once the correlation function has been selected and the best $θ_{k}$ estimated, the Kriging model can be used to predict the response $y (X)$ at an untried location $x$ using

\hat{y} (X) = \hat{β} + r^{T} (X) R^{- 1} (Y - f \hat{β}),

where $\hat{y}$ is the vector of estimated response values at each sample point, $f$ is the vector with values of the trend function evaluated at each sample point, $\hat{β}$ is a constant, and $r^{T} (X)$ is the vector of correlation values between the untried location $x$ and the sample data points.

The constant $\hat{β}$ can be estimated using the equation

\hat{β} = {(f^{T} R^{- 1} f)}^{- 1} f^{T} R^{- 1} Y .

The estimate of the variance is

σ^{2} = \frac{{(Y - f \hat{β})}^{T} R^{- 1} (Y - f \hat{β})}{n_{s}} .

The maximum likelihood estimate (i.e., best) for $θ_{k}$ is obtained by maximizing the likelihood estimate given by

\frac{n_{s} \ln (σ^{2}) + \ln | R |}{2}

Isight also supports creation of Isotropic approximations with the Kriging model. Isotropic approximation, as the name implies, is used when all independent variables behave similarly; Isight consequently assumes that all $θ_{k}$ values are identical. Because only one optimum theta value is searched for, isotropic fit is usually faster.

Depending on the number of input parameters, the number of design points, and the number of responses (outputs) of the Kriging model, the model building process can be time consuming. As the size of the matrices increases, the CPU power required for manipulating the matrices grows exponentially. Therefore, generating a good Kriging model that uses many design points can take a substantial amount of time even after all data points are analyzed.

The quality of a Kriging model depends on the location of the sample points in the design space. The Kriging model has been observed to perform best with space-filling designs where sample points are placed far apart. When points are clustered together, the matrices used in fitting the Kriging model become ill-conditioned, resulting in a poor fit. To avoid ill-conditioning, you can filter points from the sample based on distance. All points that are closer than a value called the Smoothing Filter are removed from the sample set before fitting. Isight uses other numerical techniques internally to improve the performance and robustness of the approximation.

Once a Kriging model has been built and deemed sufficiently accurate, it is ready to be used in design analysis. The time taken for fitting a Kriging model is observably larger than that required for other interpolation techniques such as Radial Basis Functions. However, the prediction times are comparable.