10 Surrogate-Based Optimization

A surrogate model, also known as a response surface model or metamodel, is an approximate model of a functional output that represents a “curve fit” to some underlying data. The goal of a surrogate model is to build a model that is much faster to compute than the original function, but that still retains sufficient accuracy away from known data points.

Surrogate-based optimization (SBO) performs optimization using the surrogate model, as shown in Figure 10.1. When used in optimization, the surrogate might define the full optimization model (i.e., the inputs are design variables, and the outputs are objective and constraint functions), or the surrogate could be just a component of the overall model. SBO is more targeted than the broader field of surrogate modeling. Instead of aiming for a globally accurate surrogate, SBO just needs the surrogate model to be accurate enough to lead the optimizer to the true optimum.

In SBO, the surrogate model is usually improved during optimization as needed but can sometimes be constructed beforehand and remain fixed during optimization. Some optimization algorithms interrogate both the surrogate model and the original model, an approach that is sometimes called surrogate-assisted optimization.

10.1 When to Use a Surrogate Model¶

There are various scenarios for which surrogate models are helpful. One scenario is when the original model is computationally expensive. Surrogate models can be queried with minimal computational cost, but constructing them requires multiple evaluations of the original model. Suppose the number of evaluations needed to build a sufficiently accurate surrogate model is less than that needed to optimize the original model directly. In that case, SBO may be a worthwhile option. Constructing a surrogate model becomes even more compelling when it is reused in multiple optimizations.

Surrogate modeling can be effective in handling noisy models because they create a smooth representation of noisy data. This can be particularly advantageous when using gradient-based optimization.

One scenario that leads to both expensive evaluation and noisy output is experimental data. When the model data are experimental and the optimizer cannot query the experiment in an automated way, we can construct a surrogate model based on the experimental data. Then, the optimizer can query the surrogate model in the optimization.

Surrogate models are also helpful when we want to understand the design space, that is, how the objective and constraints (outputs) vary with respect to the design variables (inputs). By constructing a continuous model over discrete data, we obtain functional relationships that can be visualized more effectively.

When multiple sources of data are available, surrogate models can fuse the data to build a single model. The data could come from numerical models with different levels of fidelity or experimental data. For example, surrogate models can calibrate numerical model data using experimental data. This is helpful because experimental data is usually much more scarce than numerical data. The same reasoning applies to low- versus high-fidelity numerical data.

One potential issue with surrogate models is the curse of dimensionality, which refers to poor scalability with the number of inputs. The larger the number of inputs, the more model evaluations are needed to construct a surrogate model that is accurate enough. Therefore, the reasons for using surrogate models cited earlier might not be enough if the optimization problem has a large number of design variables.

The SBO process is shown in Figure 10.2. First, we use sampling methods to choose the initial points to evaluate the function or conduct experiments. These points are sometimes referred to as training data. Next, we build a surrogate model from the sampled points. We can then perform optimization by querying the surrogate model. Based on the results of the optimization, we include additional points in the sample and reconstruct the surrogate (infill).

We repeat this process until some convergence criterion or a maximum number of iterations is reached. In some procedures, infill is omitted; the surrogate is entirely constructed upfront and not subsequently updated.

The optimization step can be performed using any of the methods we covered previously. Because surrogate models are smooth and their gradients are easily computed, gradient-based optimization is preferred (see Chapter 4). However, some surrogate models can be highly multimodal, in which case a global search is preferred, either using gradient-based with multistart (see Tip 4.8) or a global gradient-free method (see Chapter 7).

This chapter discusses sampling, constructing a surrogate, and infill with some associated optimization strategies. We devote separate sections to two surrogate modeling methods that are more involved and widely used: kriging and deep neural nets. Many of the concepts discussed in this chapter have a wide range of applications beyond optimization.

10.2 Sampling¶

Sampling methods, also known as sampling plans, select the evaluation points to construct the initial surrogate. These evaluation points must be chosen carefully. A straightforward approach is full factorial sampling, where we discretize each dimension and evaluate at all combinations of the resulting grid. This is not efficient because it scales exponentially with the number of input variables.

Example 10.1 highlights one of the significant challenges of sampling methods: the curse of dimensionality.

For SBO, even with better sampling plans, using a large number of variables is costly. We need to identify the most important or most influential variables. Knowledge of the particular domain is helpful, as is exploring the magnitude of the entries in a gradient vector (Chapter 6) across multiple points in the domain. We can use various strategies to help us decide which variables matter most, but for our purposes, we assume that the most influential variables have already been determined so that the dimensionality is reasonable. Having selected a set of variables, we are now interested in sampling methods that characterize the design space of interest more efficiently than full factorial sampling.

In addition to their use in SBO, the sampling strategies discussed in this section are useful in many other applications, including various applications discussed in this book: initializing a genetic algorithm (Section 7.6), particle swarm optimization (Section 7.7) or a multistart gradient-based algorithm (Tip 4.8), or choosing the points to run in a Monte Carlo simulation (Section 12.3.3). Because the function behavior at each sample is independent, we can efficiently parallelize the evaluation of the functions.

10.2.1 Latin Hypercube Sampling¶

Latin hypercube sampling (LHS) is a popular sampling method that is built on a random process but is more effective and efficient than pure random sampling. Random sampling scales better than full factorial searches, but it tends to exhibit clustering and requires many points to reach the desired distribution (i.e., the law of large numbers). For example, Figure 10.3 compares 50 randomly generated points across uniform distributions in two dimensions versus Latin hypercube sampling. In random sampling, each sample is independent of past samples, but in LHS, we choose all samples beforehand to ensure a well-spread distribution.

To describe the methodology, consider two random variables with bounds, whose design space we can represent as a square.

Say we wanted only eight samples; we could divide the design space into eight intervals in each dimension, generating the grid of cells shown in Figure 10.4.

A full factorial search would identify a point in each cell, but that does not scale well. To be as efficient as possible and still cover the variation, we would want each row and each column to have one sample in it. In other words, the projection of points onto each dimension should be uniform. For example, the left side of Figure 10.5 shows the projection of a uniform LHS onto each dimension. We see that the points create a uniformly spread histogram.

The concept where one and only one point exists in any given row or column is called a Latin square, and the generalization to higher dimensions is called a Latin hypercube. There are many ways to achieve this, and some choices are better than others. Consider the sampling plan shown on the left of Figure 10.6. This plan meets our criteria but clearly does not fill the space and likely will not capture the relationships between design parameters well. Alternatively, the right side of Figure 10.6 has a sample in each row and column while also spanning the space much more effectively.

LHS can be posed as an optimization problem where we seek to maximize the distance between the samples.

The constraint is that the projection on each axis must follow a chosen probability distribution. The specified distribution is often uniform, as in the previous examples, but it could also be any distribution, such as a normal distribution, as shown on the right side of Figure 10.5. This optimization problem does not have a unique solution, so random processes are used to determine the combination of points. Additionally, points are not usually placed in cell centers but at a random location within a given cell to allow for the possibility of reaching any point in the domain. The advantage of the LHS approach is that rather than relying on the law of large numbers to fill out our chosen probability distributions, we enforce it as a constraint. This method may still require many samples to characterize the design space accurately, but it usually requires far fewer than pure random sampling.

Instead of defining LHS as an optimization problem, a much simpler approach is typically used in which we ensure one sample per interval, but we rely on randomness to choose point combinations. Although this does not necessarily yield a maximum spread, it works well in practice and is simple to implement. Before discussing the algorithm, we discuss how to generate other distributions besides just uniform distributions.

We can convert from uniformly sampled points to an arbitrary distribution using a technique called inversion sampling.

Assume that we want to generate samples $x$ from an arbitrary probability density function (PDF) $p(x)$ or, equivalently, from the corresponding cumulative distribution function (CDF) $P(x)$.^[1]The probability integral transform states that for any continuous CDF, $y = P(x)$, the variable $y$ is uniformly distributed (a simple proof, but it is not shown here to avoid introducing additional notation). The procedure is to randomly sample from a uniform distribution (e.g., generate $y$), then compute the corresponding $x$ such that $P(x) = y$, which we denote as $x = P^{-1} y$. This latter step is known as an inverse CDF, a percent-point function, or a quantile function. This process is depicted in Figure 10.7 for a normal distribution. This same procedure allows us to use LHS with any distribution, simply by generating the samples on a uniform distribution.

A typical algorithm is described in Algorithm 10.1. For each axis, we partition the CDF in $n_s$ evenly spaced regions (evenly spaced along the CDF, which means that each region is equiprobable). We generate a random number within each evenly spaced interval, where 0 corresponds to the bottom of the interval and 1 to the top. We then evaluate the inverse CDF as described previously so that the points match our specified distribution (the CDF for a uniform distribution is just a line $P(x) = x$, so the output is not changed). Next, the column of points for that axis is randomly permuted. This process is repeated for each axis according to its specified probability distribution.

An example using Algorithm 10.1 for eight points is shown in Figure 10.8. In this example, we use a uniform distribution for $x_1$ and a normal distribution for $x_2$. There is one point in each equiprobable interval. As stated before, randomness does not necessarily ensure a good spread, but optimizing the spread is difficult because the function is highly multimodal. Instead, to encourage high spread, we could generate multiple Latin hypercube samples with this algorithm and select the one with the largest sum of the distance between points.

10.2.2 Low-Discrepancy Sequences¶

Low-discrepancy sequences generate deterministic sequences of points that are well spread. Each new point added in the sequence maintains low discrepancy—discrepancy refers to the variation in point density throughout the domain. Hence, a low-discrepancy set is close to even density (i.e., well spread). These sequences are called quasi-random because they often serve as suitable replacements for applications that use random sequences, but they are not random or even pseudo-random.

An advantage of low-discrepancy sequences over LHS is that most of the approaches do not require selecting all the samples beforehand. These methods generate deterministic sequences; in other words, we generate the same sequence of points whether we choose them beforehand or add more later. This property is particularly advantageous in iterative procedures. We may choose an initial sampling plan and add more well-spread points to the sample later. This is not necessarily an advantage for the methods of this chapter because the optimization drives the selection of new points rather than continuing to seek spread-out samples. However, this feature is useful for other applications, such as quadrature, Monte Carlo simulations, and other problems where an iterative sampling process is used to determine statistical convergence (see Section 12.3). Low-discrepancy sequences add more points that are well spread without having to throw out the existing samples. Even in non-iterative procedures, these sampling strategies can be a useful alternative.

Several of these sequences are built on generalizations of the one-dimensional van der Corput sequence to more than one dimension. Such sequences are defined by representing an integer $i$ in a given integer base $b$ (the van der Corput sequence is always base 2):

If the base is 2, this is just a standard binary sequence (Section 7.6.1). After determining the relevant coefficients ($a_j$), the $i\text{th}$ element of the sequence is

An algorithm to generate an element in this sequence, also known as a radical inverse function for base $b$, is given in Algorithm 10.2.

For base 2, the sequence is as follows:

The interval is divided in half, and then each subinterval is also halved, with new points spreading out across the domain (see Figure 10.9).

Similarly, for base 3, the interval is split into thirds, then each subinterval is split into thirds, and so on:

Halton Sequence¶

A Halton sequence uses pairwise prime numbers (larger than 1) for the base of each dimension of the problem.^[2]The $i\text{th}$ point in the Halton sequence is

where the $b_j$ set is pairwise prime. As an example in two dimensions, Figure 10.10 shows 30 generated points of the Halton sequence where $x_1$ uses base 2, and $x_2$ uses base 3, and then a subsequent 20 generated points are added (in another color), showing the reuse of existing points.

If the dimensionality of the problem is high, then some of the base combinations lead to points that are highly correlated and thus undesirable for a sampling plan. For example, the left of Halton sequence with base 17 for $x_1$ and base 19 for $x_2$. shows 50 generated points where $x_1$ uses base 17, and $x_2$ uses base 19. To avoid this issue, we can use a scrambled Halton sequence.

Scrambling can be accomplished by generating a permutation array containing a random permutation of the integers $p = [0, 1, \ldots, b-1]$. Then, rather than using the integers $a$ directly in Equation 10.2, we use the entries of $a$ as the indices of the permutation array. If $p$ is the permutation array, we have:

The permutation array is fixed for all digits $a$ and for all $n_p$ points in the domain. The right side of Halton sequence with base 17 for $x_1$ and base 19 for $x_2$. shows the same example (with base 17 and base 19) but with scrambling to weaken the strong correlations.

Hammersley Sequence¶

The Hammersley sequence is closely related to the Halton sequence. However, it provides better spacing if we know beforehand the number of points ($n_p$) that we are going to use. This approach only needs $n_x - 1$ bases (still pairwise prime) because the first dimension uses regular spacing:

Because this sequence needs to know the number of points ($n_p$) beforehand, it is less useful for iterative procedures.

However, the implementation is straightforward, and so it may still be a useful alternative to LHS. Figure 10.12 shows 50 points generated from a Hammersley sequence where the $x_2$-axis uses base 2.

Other Sequences¶

A wide variety of other low-discrepancy sequences exist. The Faure sequence is similar to the Halton, but it uses the same base for all dimensions and uses permutation scrambling for each dimension instead.12 Sobol sequences use base 2 sequences but with a reordering based on “direction numbers”.3 Niederreiter sequences are effectively a generalization of Sobol sequences to other bases.4

10.3 Constructing a Surrogate¶

Once sampling is completed, we have a list of data points, often called training data:

where $x^{(i)}$ is an input vector from the sampling plan, and $f^{(i)}$ contains the corresponding outputs from evaluating the model: $f^{(i)} = f\left(x^{(i)}\right)$. We seek to construct a surrogate model from this data set. Surrogate models can be based on physics, mathematics, or a combination of the two.

Incorporating known physics into a model is often desirable to improve model accuracy. However, functional relationships are unknown for many complex problems, and a data-driven mathematical model can be more effective.

Surrogate-based models can be based on interpolation or regression, as illustrated in Figure 10.13. Interpolation builds a function that exactly matches the provided training data. Regression models do not try to match the training data points exactly; instead, they minimize the error between a smooth trend function and the training data. The nature of the training data can help decide between these two types of surrogate models. Regression is particularly useful when the data are noisy. Interpolatory models may produce undesirable oscillations when fitting the noise. In contrast, regression models can find a smooth function that is less sensitive to the noise. Interpolation is useful when the data are highly multimodal (and not noisy). This is because a regression model may smooth over variations that are actually physical, whereas an interpolatory model can accurately capture those variations.

There are two main steps involved in either type of surrogate model. First, we select a set of basis functions, which represent the form for the model.

Second, we determine the model parameters that provide the best fit to the provided data. Determining the model parameters is an optimization problem, which we discuss first. We discuss linear regression and nonlinear regression, which are techniques for choosing model parameters for a given set of basis functions. Next, we discuss cross validation, which is a critical technique for selecting an appropriate model form. Finally, we discuss common basis functions.

10.3.1 Linear Least Squares Regression¶

A linear regression model does not mean that the surrogate is linear in the input variables but rather that the model is linear in its coefficients (i.e., linear in the parameters we are estimating). For example, the following equation is a two-dimensional linear regression model, where we use $\hat{f}$ to represent our estimated model of the function $f$:

This function is highly nonlinear, but it is classified as a linear regression model because the regression seeks to choose the appropriate values for the coefficients $w_i$ (and the function is linear in $w$).

A general linear regression model can be expressed as

where $w$ is a vector of weights, and $\psi$ is a vector of basis functions. In this section, we assume that the basis functions are provided. In general, the basis functions can be any set of functions that we choose (and typically they are nonlinear). It is usually desirable for these functions to be orthogonal.

The coefficients are chosen to minimize the error between our predicted function values $\hat{f}$ and the actual function values $f^{(i)}$. Because we want to minimize both positive and negative errors, we minimize the sum of the square of the errors (or a weighted sum of squared errors):^[3]

The solution to this optimization problem is called a least squares solution. If the regression model is linear, we can simplify this objective and solve the problem analytically. Recall that $\hat{f} = \psi^\intercal w$, so the objective can be written as

We can express this in matrix form by defining the following:

The matrix $\Psi$ is of size $(n_s \times n_w)$, where $n_s$ is the number of samples, $n_w$ the number of parameters in $w$, and $n_s \ge n_w$. This means that there should be more equations than unknowns or that we have sampled more points than the number of coefficients we need to estimate. This should make sense because our surrogate function is only an assumed form and generally not an exact fit to the actual underlying function. Thus, we need more data to create a good fit.

Then the optimization problem can be written in matrix form as:

Expanding the squared norm (i.e., $\|x\|_2^2 = x^\intercal x$) gives

We can omit the last term from the objective because our optimization variables are $w$, and the last term has no $w$ dependence:

This fits the general form for an unconstrained quadratic programming (QP) problem, as shown in Section 5.5.1:

where

Recall that an equality constrained QP (of which unconstrained is a subset) has an analytic solution as long as the QP is positive definite. In our case, we can show that $Q$ is positive definite as long as $\Psi$ is full rank:

This is not surprising because the objective is a sum of squared values. Referring back to the solution in Section 5.5.1, and removing the portions associated with the constraints, the solution is

In our case, this becomes

After simplifying, we have an analytic solution for the weights:

We sometimes express the linear relationship in Equation 10.12 as $\Psi w = f$, although the case where there are more equations than unknowns does not typically have a solution (the problem is overdetermined). Instead, we seek the solution that minimizes the error $\|\Psi w - f\|^2$, that is, Equation 10.23. The quantity $\Psi^\dagger = \left(\Psi^\intercal \Psi \right)^{-1} \Psi^\intercal$ is called the pseudoinverse of $\Psi$ (or more specifically, the Moore–Penrose pseudoinverse), and thus we can write Equation 10.23 in the more compact form

This allows for a similar form to solving a linear system of equations where an inverse would be used instead. In solving both a linear system and the linear least-squares equation (Equation 10.23), we do not explicitly invert a matrix. For linear least squares, a QR factorization is commonly used for improved numerical conditioning as compared to solving Equation 10.23 directly.

A common variation of this approach is to use regularized least squares, which adds a term in the objective. The new objective is

where $\mu$ is a weight assigned to the second term. This second term attempts to reduce the magnitudes of the entries in $w$ while balancing the fit in the first term. This approach is particularly beneficial if the data contain strong outliers or are particularly noisy. The rationale for this approach is that we may want to accept a higher error (quantified by the first term) in exchange for smaller values for the coefficients. This generally leads to simpler, more generalizable models (e.g., by reducing the influence of some terms). A related extension uses a second term of the form $\| w - w_0 \|_2^2$. The idea is that we want a good fit, while maintaining parameters that are close to some known nominal values $w_0$.

A regularized least squares problem can be solved with the same linear least squares approach. We can write the previous problem using concatenated matrices and vectors:

This is of the same linear form as before ($\| Aw - b \|^2$), so we can reuse the solution (Equation 10.23):

For linear least squares (with or without regularization), we have seen that the optimization problem of determining the appropriate coefficients can be found analytically. We can also add linear constraints to the problem (equality or inequality), and the optimization remains a QP. In that case, the problem is still convex. Although it does not generally have an analytic solution, we can still quickly find the global optimum. This topic is discussed in Section 11.3.

10.3.2 Maximum Likelihood Interpretation¶

This section presents an alternative motivation for the sum of squared error approach used in the previous section. It is somewhat of a diversion from the present discussion, but it will be helpful in several results later in this chapter. In the previous section, we assumed linear models of the form

where we use $x$ for simplicity in writing instead of $\psi(x^{(i)})$. The derivation remains the same for any arbitrary function of $x$. The function $\hat{f}$ is just a model, so we could say that it is equal to our actual observations $f^{(i)}$ plus an error term:

where $\varepsilon$ captures the error associated with the $i{\text{th}}$ data point. We assume that the error is normally distributed with mean zero and a standard deviation of $\sigma$:

The use of Gaussian uncertainty can be motivated by the central limit theorem, which states that for large sample sizes, the sum of random variables tends toward a Gaussian distribution regardless of the original distribution associated with each variable. In other words, the sample distribution of the sample mean is approximately Gaussian. Because we assume the error terms to be the sum of various random, independent perturbations, then by the central limit theorem, we expect the errors to be normally distributed.

We now substitute Equation 10.29 into Equation 10.30 to show the probability of $f$ conditioned on $x$ and parameterized by $w$:

Once we include all the data points $x^{(i)}$, we would like to compute the probability of observing $f$ conditioned on the inputs $x$ for a given set of parameters in $w$. We call this the likelihood function $L(w)$. In this case, assuming all errors are independent, the total probability for observing the outputs is the product of the probability of observing each output:

Now we can pose this as an optimization problem where we wish to find the parameters $w$ that maximize the likelihood function; in other words, we maximize the probability that our model is consistent with the observed data. Because the objective is a product of multiple terms, it is helpful to take the logarithm of the objective. Maximizing $L$ or maximizing $\ell = \ln(L)$ does not change the solution to the problem but makes it easier to solve. We call this the log likelihood function:

The first term has no dependence on $w$, and so when optimizing $\ell(w)$; it is just a scalar term that can be removed as follows:

Similarly, the denominator of the second term has no dependence on $w$ and is just a scalar that can also be removed:

Thus, maximizing the log likelihood function (maximizing the probability of observing the data) is equivalent to minimizing the sum of squared errors (the least-squares formulation). This derivation provides another motivation for using the sum of squared errors in regression.

10.3.3 Nonlinear Least Squares Regression¶

A surrogate model can be nonlinear in the coefficients. For example, building on the simple function shown earlier in Equation 10.9 we can add the coefficients $w_6$ and $w_7$ as follows:

The addition of $w_6$ and $w_7$ makes this function nonlinear in the coefficients. We can still estimate these parameters, but not analytically.

Equation 10.11 is still relevant; we still seek to minimize the sum of the squared errors:

For general nonlinear regression models, we cannot write a more specific form for $\hat{f}$ as we could for the linear case, so we leave the objective as it is.

This is a nonlinear least-squares problem. The optimization problem is unconstrained, so any of the methods from Chapter 4 apply. We could also easily add constraints, for example, bounds on parameters, known relationships between parameters, and so forth, and use the methods from Chapter 5.

In contrast to the linear case, we need to provide a starting point, our best guess for the parameters, and we may need to deal with scaling, noise, multimodality, or any of the other potential challenges of general nonlinear optimization. Still, this is a relatively straightforward problem within the broader realm of engineering optimization problems.

Although the methods of Chapter 4 can be used if the problem remains unconstrained, there are more specialized methods available that take advantage of the specific structure of the problem. One popular approach to solving the nonlinear least-squares problem is the Levenberg–Marquardt algorithm, which we discuss in this section.

As a stepping stone towards the Levenberg–Marquardt algorithm, we first derive the Gauss–Newton algorithm, which is a modification of Newton’s method (Section 3.8) for solving nonlinear least-squares problems. One way to think of this algorithm is as an iterative linearization of the residual. Once it is linearized, we can apply the same methods we derived for linear least squares. We linearize the residual $r = \hat{f}(w) - f$ at iteration $k$ as

where $\Delta w$ is the step and the Jacobian is

After the linearization, the objective becomes

This is now the same form as linear least squares (Equation 10.14), so we can reuse its solution (Equation 10.23) to solve for the step

We now have an update formula for the coefficients at each iteration:

An alternative derivation for this formula is to consider a Newton step for an unconstrained optimizer. The objective is $e = \sum_i r_i^2$, and the formula for a Newton step (Section 4.4.3) is

The gradient is

or in matrix form:

The Hessian in index notation is

We can write it in matrix form as follows:

If we neglect the second term in the Hessian, then the Newton update is:

which is the same update as before.

Thus, another interpretation of this method is that a Gauss–Newton step is a modified Newton step where the second derivatives of the residual are neglected (and thus, a quasi-Newton approach to estimate second derivatives is not needed). This method is particularly effective near convergence because as $r \rightarrow 0$ (i.e., as we approach the solution to our residual minimization), the neglected term also approaches zero. The appeal of this approach is that we can often obtain an accurate prediction for the Hessian using only the first derivatives because of the known structure of the objective.

When the second term is not small, then the Gauss–Newton step may be too inaccurate. We could use a line search, but the Levenberg–Marquardt algorithm utilizes a different strategy. The idea is to regularize the problem as discussed in the previous section or, in other words, provide the ability to dampen the steps as needed. Each linearized subproblem becomes

Recall that the solution to this problem (see Equation 10.27) is

If $\mu = 0$, then we retain the Gauss–Newton step. Conversely, as $\mu$ becomes large, so that the $J_r^\intercal J_r$ is negligible, the step becomes

This is precisely the steepest-descent direction for our objective (see Equation 10.49), although with a small magnitude because $\mu$ is large. The parameter $\mu$ provides some control for directions ranging between Gauss–Newton and steepest descent.

The Levenberg–Marquardt algorithm has been revised to improve the scaling for components of the gradient that are small. The second minimization term weights all parameters equally. The scaling can be improved by multiplying by a diagonal matrix in the regularization as follows:

where $D$ is defined as

This matrix scales the objective by the diagonal elements of the Hessian. Thus, when $\mu$ is large, and the direction tends toward the steepest descent, the components of the gradient are scaled by the curvature, which reduces the amount of zigzagging.

The solution to the minimization problem of Equation 10.56 is

Finally, we describe one of the possible heuristics for selecting and updating the damping parameter $\mu$. After a successful step (a sufficient reduction in the objective), $\mu$ is increased by a factor of $\rho$. Conversely, an unsuccessful step is rejected, and $\mu$ is reduced by a factor ($\mu = \mu / \rho$).

Rather than returning a scalar objective ($\sum_i r_i^2$), the user function should return a vector of the residuals because that vector is needed in the update steps (Equation 10.58). A potential convergence metric is a tolerance on objective value changes between subsequent iterations. The full procedure is described in Algorithm 10.3.

10.3.4 Cross Validation¶

The other important consideration for developing a surrogate model is the choice of the basis functions in $\psi$. In some instances, we may know something about the model behavior and thus what type of basis functions should be used, but generally, the best way to determine the basis functions is through cross validation. Cross validation is also helpful in characterizing error, even if we already have a chosen set of basis functions. One of the reasons we use cross validation is to prevent overfitting. Overfitting occurs when we have too many degrees of freedom and closely fit a given set of data, but the resulting model has a poor predictive ability. In other words, we are fitting noise. The following example illustrates this idea with a one-dimensional function.

The solution to the overfitting problem highlighted in Example 10.5 is cross validation. Cross validation means that we use one set of data for training (creating the model) and a different set of data for assessing its predictive error. There are many different ways to perform cross validation; we describe two. Simple cross validation is illustrated in Figure 10.17 and consists of the following steps:

1. Randomly split your data into a training set and a validation set (e.g., a 70–30 split).

2. Train each candidate model (the different options for $\psi$) using only the training set, but evaluate the error with the validation set. The error on previously unseen data is called the generalization error ($e_g$ in Figure 10.17).

3. Choose the model with the lowest generalization error, and optionally retrain that model using all of the data.

An alternative option that is more involved but uses the data more effectively is called $k$-fold cross validation. It is particularly advantageous when we have a small data set where we cannot afford to leave much out. This procedure is illustrated in Figure 10.18 and consists of the following steps:

1. Randomly split your data into $n$ sets (e.g., $n$ = 10).

2. Train each candidate model using the data from all sets except one (e.g., 9 of the 10 sets) and use the remaining set for validation. Repeat for all $n$ possible validation sets and average the performance.

3. Choose the model with the lowest average generalization error. Optionally, retrain with all the data.

The extreme version of this process, when training data are very limited, is leave-one-out cross validation (i.e., each testing subset consists of one data point).

10.3.5 Common Basis Functions¶

Although cross validation can help us find the lowest generalization error among a provided set of basis functions, we still need to determine what sets of options to consider. This selection is crucial because our model is only as good as the available options, but increasing the number of options increases computational time. The possibilities for basis functions are as numerous as the types of function. As stated before, it is generally desirable that they form an orthogonal set. We focus on a few commonly used functions.

Polynomials¶

Polynomials, of which we have already seen a few examples, are useful in many applications. However, we typically use low-order polynomials for regression because high-order polynomials rarely generalize well. Polynomials can be particularly effective in cases where a knowledge of the physics suggests them to be an appropriate choice (e.g., drag varies quadratically with speed)

Because a lot of structure is already built into the model form, fewer data points are needed to create a reasonable model (e.g., a quadratic function in $n$ dimensions needs at least $n(n+1)/2 + n + 1$ points, so this amounts to 6 points in two dimensions, 10 points in three dimensions, and so on).

Radial Basis Functions¶

Another common type of basis function is a radial basis function. Radial basis functions are functions that depend on the distance from some center point and can be written as follows:

where $c$ is the center point, and $r$ is the radius about the center point. Although the center points can be placed anywhere, we usually choose the sampling data as centering points:

This is often a useful choice because it captures the idea that our ability to predict function behavior is related to how close we are to known function values (in other words, nearby points are more highly correlated). This form naturally lends itself to interpolation, although regularization can be added to allow for regression. Polynomials are often combined with radial basis functions because the polynomial can capture global function behavior, while the radial basis functions can introduce modifications to capture local behavior.

One popular radial basis function is the Gaussian basis:

where $\theta_j$ are the model parameters. One of the forms of kriging discussed in the following section can be viewed as a radial basis function model with a Gaussian basis.

10.4 Kriging¶

Kriging is a popular surrogate modeling technique that can build approximations of highly nonlinear engineering simulations. We may not have a simple parametric form for such simulations that we can use with regression and expect a good fit.

Instead of tuning the parameters of a functional form that describes what the function is, kriging tunes the parameters of a statistical model that describes how the function behaves.

The kriging statistical model that approximates $f$ consists of two terms: a function $\mu(x)$ that is meant to capture some of the function behavior and a random variable $Z(x)$. Thus, we can write the kriging model as

When we evaluate the function we want to approximate at point $x$, we get a scalar value $f(x)$. In contrast, when we evaluate the stochastic process (Equation 10.62) at $x$, we get a random variable $F(x)$ that has a normal distribution with mean $\mu$ and variance $\sigma^2$. Although we wrote $\mu$ as a function of $x$, most kriging models consider this to be constant because the random variable term alone is effective in capturing the function behavior.

For the rest of this section, we discuss the case with constant $\mu$, which is called ordinary kriging. Kriging is also referred to as Gaussian process interpolation, or more generally in the regression case discussed later in this section, as Gaussian process regression.

The power of the statistical model lies in how it treats the correlation between the random variables. Although we do not know the exact form of the error term $Z(x)$, we can still make some reasonable assumptions about it. Consider two points in a sampling plan, $x^{(i)}$ and $x^{(j)}$, and the corresponding terms, $Z\left(x^{(i)}\right)$ and $Z\left(x^{(j)}\right)$. Intuitively, we expect $Z\left(x^{(i)}\right)$ to be close to $Z\left(x^{(j)}\right)$ whenever $x^{(i)}$ is close to $x^{(j)}$. Therefore, it seems reasonable to assume that the correlation between $Z\left(x^{(i)}\right)$ and $Z\left(x^{(j)}\right)$ is a function of the distance between the two points.

In kriging, we assume that this correlation is given by a kernel function $K\left( x^{(i)},\, x^{(j)}\right)$:

As a matrix, the kernel is represented as $K_{ij} = K\left(x^{(i)}, x^{(j)}\right)$. Various kernel functions are used with kriging.^[5]The most commonly used kernel function is

where $\theta_l \ge 0$, $0 \le p_l \le 2$, and $n_d$ is the number of dimensions (i.e., the length of the vector $x$). If every $p_l = 2$, this becomes a Gaussian kernel.

Let us examine how the statistical model $F$ defined in Equation 10.62 captures the typical behavior of the function $f$.

The parameter $\mu$ captures the typical value, and $\sigma^2$ captures the expected variance. The kernel (or correlation) function (Equation 10.64) implicitly models continuous functions. If $f$ is continuous, we know that, as $|x^{(i)}-x^{(j)}| \rightarrow 0$, then $|f(x^{(i)})-f(x^{(j)})| \rightarrow 0$. This is captured in the kernel function because as $|x^{(i)}-x^{(j)}| \rightarrow 0$, the correlation approaches 1. The parameter $\theta_l$ captures how active the function $f$ is in the $l$th coordinate direction. A unit difference in variable $l$ ($|x^{(i)}_l - x^{(j)}_l| = 1$) has a more significant impact on the correlation when $\theta_l$ is large. The exponent $p_l$ describes the smoothness of the function in the $l$th coordinate direction. Values of $p_l$ close to 2 produce smooth functions, whereas values closer to zero produce functions with more variation.

Kriging surrogate modeling involves two main steps.

The first step consists of using the data to estimate the statistical model parameters $\mu$, $\sigma^2$, $\theta_1,\ldots,\theta_{n_d}$, and $p_1,\ldots,p_{n_d}$. The second step consists of making predictions using the statistical model and these estimated parameter values.

The parameter estimation uses the same maximum likelihood approach from Section 10.3.2, but now it is more complicated. Let us denote the random variable as $F^{(i)} \equiv F\left(x^{(i)}\right)$ and the vector of random variables as $F = \left[ F^{(1)}, \ldots, F^{(n_s)}\right]$, where $n_s$ is the number of samples. Similarly, $f^{(i)} \equiv f\left(x^{(i)}\right)$ and the vector of observed function values is $f \equiv [f^{(1)},\dots,f^{(n_s)}]$. Using this notation, we can say that the vector $F$ is jointly normally distributed. This is also known as a multivariate Gaussian distribution.^[6]The probability density function (PDF) (the likelihood that ${F}={f}$) is

where $e$ is a vector of 1s with size $n_s$, and $|\Sigma|$ is the determinant of the covariance,

The covariance between two elements $F^{(i)}$ and $F^{(j)}$ of ${F}$ is related to correlation by the following definition:^[7]

We assume stationarity of the second moment, that is, the variance $\sigma^2$ is constant in the domain.

The statistical model parameters $\theta_1,\ldots,\theta_{n_d}$ and $p_1,\ldots,p_{n_d}$ enter the likelihood (Equation 10.65) via their effect on the kernel $K$ (Equation 10.64) and hence on the covariance matrix $\Sigma$ (Equation 10.66).

We estimate the parameters using the maximum log likelihood approach from Section 10.3.2; that is, we maximize the probability of observing our data ${f}$ conditioned on the parameters $\mu$ and $\Sigma$. Using a PDF (Equation 10.65) where $\mu$ is constant and the covariance is $\Sigma = \sigma^2 K$, yields the following likelihood function:

We now need to find the parameters $\mu$, $\sigma$, $\theta_i$, and $p_i$ that maximize this likelihood function, that is, maximize the probability of our observations ${f}$.

As before, we take the logarithm to form the log likelihood function:

We can maximize part of this term analytically by taking derivatives with respect to $\mu$ and $\sigma$, setting them equal to zero, and solving for their optimal values to obtain:

We now substitute these values back into the log likelihood function (Equation 10.68), which yields

This function, also called the concentrated likelihood function, only depends on the kernel $K$, which depends on $\theta$ and $p$.

We cannot solve for optimal values of $\theta$ and $p$ analytically. Instead, we rely on numerical optimization to maximize Equation 10.71. Because $\theta$ can vary across a broad range, it is often better to search using logarithmic scaling. Once we solve that optimization problem, we compute the mean and variance in Equation 10.69 and .

Now that we have a fitted model, we can make predictions at new points where we have not sampled. We do this by substituting $x_p$ into a formula called the kriging predictor. The formula is unique, but there are many ways to derive it.

One way to derive it is to find the function value at $x_p$ that is the most consistent with the behavior of the function captured by the fitted kriging model.

Let $f_p$ be our guess for the value of the function at $x_p$. One way to assess the consistency of our guess is to add $(x_p, f_p)$ as an artificial point to our training data (so that we now have $n_s+1$ points) and estimate the likelihood using the parameters from our fitted kriging model. The likelihood of this augmented data can now be thought of as a function of $f_p$: high values correspond to guessed values of $f_p$ that are consistent with function behavior captured by the fitted kriging model. Therefore, the value of $f_p$ that maximizes the likelihood of this augmented data set is a natural way to predict the value of the function. This is an optimization problem with a closed-form solution, and the corresponding formula is the kriging predictor.

Now we outline the derivation of the kriging predictor.^[8]

With the augmented point, our function values are $\bar{f} = [{f}, f_p]$, where $f$ is the $n_s$-vector of function values from the original training data. Then, the correlation matrix with the additional data point is

where $k$ is the correlation of the new point with the training data given by

The 1 in the bottom right of the augmented correlation matrix (Equation 10.72) is because the correlation of the new variable $F(x_p)$ with itself is 1. The log likelihood function with these new augmented vectors and the previously determined parameters is as follows (see Equation 10.68):

We want to maximize this function with respect to $f_p$. Because only the last term depends on $f_p$ (it is a part of $\bar{f}$) we can omit the other terms and formulate the following:

This problem can be solved analytically, yielding the mean value of the kriging prediction,

The mean square error of the kriging prediction (that is, the expected squared value of the error) is given by^[9]

One attractive feature of kriging models is that they are interpolatory and thus match the training data exactly. To see how this is true, if $x_p$ is the same as one of the training data points, $x^{(i)}$, then $k$ is just $i$th column of $K$. Hence, $K^{-1}k$ is a vector $e_i$, with all zeros except for 1 in the $i$th element. In the prediction (Equation 10.75), $k^\intercal K^{-1}=e_i^\intercal$ and so the last term is $f^{(i)}-\mu^*$, which means that $f_p = f^{(i)}$.

In the mean square error (Equation 10.76), $k^\intercal K^{-1} k$ is the same as $k^\intercal e_i$. This is the $i$th element of $k$, which is 1. Therefore, the first two terms in the brackets in Equation 10.76 cancel, and the last term is zero, yielding $\sigma^2_p=0$. This is expected; if we already sampled the point, the uncertainty about its function value should be zero.

When describing a fitted kriging model, we often refer to the standard error as the square root of this quantity (i.e., $\sqrt{\sigma_{p}^2}$). The standard error is directly related to the confidence interval (e.g., $\pm 1$ standard error corresponds to a 68 percent confidence interval).

If we can provide the gradients of the function at the training data points (in addition to the function values), we can use that information to build a more accurate kriging model. This approach is called gradient-enhanced kriging (GEK). The methodology is the same as before, except we add more observed outputs (i.e., in addition to the function values at the sampled points, we add their gradients). In addition to considering the correlation between the function values at different sampled points, the kernel matrix $K$ needs to be expanded to consider correlations between function values and gradients, gradients and function values, and among gradient components.

We can use still use equation (Equation 10.75) for the GEK predictor and equation (Equation 10.76) for the mean square error if we plug in “expanded versions” of the outputs ${f}$, the vector $k$, the matrix $K$, and the vector of 1s, $e$.

We expand the output vector to include not just the function values at the sampled points but also their gradients:

This vector is of length $n_s + n_s n_d$, where $n_d$ is the dimension of $x$. The gradients are usually provided at the same $x$ locations as the function samples, but that is not required.

Recall that the term $e\mu^*$ in Equation 10.75 for the kriging predictor represents the expected value of the random variables $F^{(1)},\ldots,F^{(n_s)}$. Now that we have expanded the outputs to include the gradients at the sampled points, the mean vector needs to be expanded to include the expected values of $\nabla F^{(i)}$, which are all zero.

We can still use $e\mu^*$ in the formula for the predictor if we use the following definition:

where 1 occurs for the first $n_s$ entries, and 0 for the remaining $n_s n_d$ entries.

The additional correlations (between function values and derivatives and between the derivatives themselves) are as follows:

Here, we use $l$ and $k$ to represent a component of a vector, and we use $K_{ij} \equiv K\left(x^{(i)}, x^{(j)}\right)$ as shorthand. For our particular kernel choice (Equation 10.64), these correlations become the following: