Bayesian Optimization with Gaussian Processes: Surrogate Models for Expensive CFD

The expensive CFD problem

A typical CFD simulation takes between 20 minutes and several hours. Genetic algorithms need thousands of evaluations to converge. If each evaluation costs 1 hour and your algorithm needs 1000 evaluations, your optimization takes 42 days. Not viable for real project timelines.

Bayesian optimization solves this: it builds a surrogate model of the objective function from evaluations already performed and uses that model to intelligently decide what point to evaluate next. You will find competitive designs with 50-100 evaluations instead of thousands. The key insight is that it is more efficient to spend computational effort on building a model that guides the search than to blindly evaluate thousands of designs.

This approach is fundamentally different from evolutionary algorithms. Instead of maintaining a population and letting it evolve through random variation, Bayesian optimization maintains a probabilistic model of the objective function and uses it to make optimal decisions about where to sample next. Each new evaluation is chosen to maximize the expected information gain about the location of the optimum — not to randomly explore the design space.

What is a Gaussian Process?

A Gaussian Process (GP) is a distribution over functions, not over vectors. It is completely defined by:

Mean function m(x): the central tendency of the function (typically m(x) = 0 after standardizing outputs).
Covariance kernel k(x, x'): how correlated the outputs of two points x and x' are. Encodes the smoothness assumption.

The fundamental property: any finite set of evaluated points has a joint Gaussian distribution with mean and covariance defined by the kernel. This allows making predictions with uncertainty — unlike a neural network or polynomial fit, a GP does not just give you a single predicted value but a full probability distribution over possible function values at any point.

Common kernels in engineering:

RBF / Squared Exponential: k(x, x') = sigma² * exp(-||x - x'||² / (2 * l²)). Assumes infinitely differentiable functions — often too smooth for engineering response surfaces.
Matern 5/2: assumes only 2 continuous derivatives, more realistic for CFD and FEA response surfaces. This is the default recommendation for engineering applications.

The choice of kernel encodes prior beliefs about the function's smoothness. An RBF kernel says "I believe this function is infinitely smooth" — which may cause the GP to be overconfident in regions between observations. The Matern 5/2 kernel says "I believe this function is rough — possibly with sharp transitions" — which produces more cautious predictions that better reflect the true behavior of physical simulations where small geometry changes can produce abrupt performance changes (e.g., stall, separation, buckling).

Below is a minimal GP implementation with the Matern 5/2 kernel, demonstrating prediction with uncertainty:

import numpy as np
from scipy.linalg import cholesky, solve_triangular

def matern52_kernel(X1, X2, length_scale=1.0, sigma=1.0):
    """Matern 5/2 kernel with a single length scale."""
    dist = np.sqrt(np.sum((X1[:, None] - X2[None, :]) ** 2, axis=-1))
    sqrt5 = np.sqrt(5.0)
    r = sqrt5 * dist / length_scale
    return sigma**2 * (1.0 + r + r**2 / 3.0) * np.exp(-r)

def gp_predict(X_train, y_train, X_test, kernel, noise=1e-6):
    """GP prediction with posterior mean and variance."""
    K = kernel(X_train, X_train)
    K += noise * np.eye(len(X_train))
    K_s = kernel(X_train, X_test)
    K_ss = kernel(X_test, X_test)
    
    L = cholesky(K, lower=True)
    alpha = solve_triangular(L.T, solve_triangular(L, y_train, lower=True))
    
    mu = K_s.T @ alpha
    v = solve_triangular(L, K_s, lower=True)
    sigma = np.diag(K_ss) - np.sum(v**2, axis=0)
    return mu, np.maximum(sigma, 0)

The Cholesky decomposition is the computational bottleneck of GP prediction, scaling as O(N³) where N is the number of training points. For the typical 50-200 points in a Bayesian optimization run, this is negligible (milliseconds). The prediction step itself is O(N) per test point. The noise parameter on the diagonal (often called the "nugget") serves a dual purpose: it ensures numerical stability of the Cholesky decomposition and accounts for small amounts of noise in the observations — even deterministic CFD simulations have small mesh-dependent variations.

Acquisition functions: deciding what to evaluate

The trained GP gives us for each point x a prediction with mean mu(x) and standard deviation sigma(x). But it does not tell us what point to evaluate next. The acquisition function quantifies numerically how promising a point is, balancing exploration (zones with high uncertainty) and exploitation (zones with low mean).

Expected Improvement (EI): The most widely used acquisition function. It computes the expected improvement over the best observed value f_best:

EI(x) = E[max(f_best - f(x), 0)]

Analytical formula: EI(x) = (f_best - mu(x)) * Phi(Z) + sigma(x) * phi(Z), where Z = (f_best - mu(x)) / sigma(x), Phi is the normal CDF and phi is the normal PDF. This elegant closed form means EI can be computed and optimized extremely efficiently — the gradient with respect to x is also analytic, enabling gradient-based optimization of the acquisition function itself.

The two terms in EI have an intuitive interpretation. The first term (f_best - mu) * Phi(Z) represents exploitation: it is large when the predicted mean is below the current best and the probability of improvement is high. The second term sigma * phi(Z) represents exploration: it is large when uncertainty is high and there is a non-negligible chance of significant improvement. The balance between the two emerges automatically from the GP's uncertainty estimates — no manual tuning parameter is required.

Upper Confidence Bound (UCB): Simpler to implement: UCB(x) = mu(x) - kappa * sigma(x) in minimization. The parameter kappa controls the exploration-exploitation balance explicitly. kappa = 2 corresponds roughly to a 95% confidence bound. Unlike EI, UCB requires manual tuning of kappa, and the optimal value depends on the problem. However, UCB's simplicity makes it attractive for quick implementations and for cases where the GP uncertainty estimates are well-calibrated.

Complete Bayesian optimization pipeline

Initial DoE: Generate 10-20 points with Latin Hypercube Sampling to efficiently cover the design space.
Evaluate CFD at all initial points. Add to dataset D.
Train GP on all data in D. Adjust kernel hyperparameters (length scales, sigma) via maximum marginal likelihood.
Optimize acquisition function: Find the point x_next that maximizes EI(x) (or minimizes UCB). Use L-BFGS or multistart random search.
Evaluate CFD at x_next and add to dataset D.
Repeat steps 3-5 until the evaluation budget is exhausted or convergence is reached.
Return the best evaluated point (not the GP minimum — always validate with real simulation).

Step 7 deserves emphasis: the GP is an approximation. Its predicted minimum may differ from the true minimum, especially in regions far from any evaluated point. Always return the best point that was actually evaluated with CFD, not the GP's predicted optimum. The GP's role is to guide the search, not to replace the simulation entirely.

The Latin Hypercube Sampling in step 1 ensures good coverage of the design space with a small number of points. Unlike simple random sampling, LHS stratifies each dimension, guaranteeing that no region of the design space is left unexplored in the initial dataset. This is critical because the GP's uncertainty estimates are only reliable in regions near training points — if large gaps exist in the initial sampling, the GP may be overly optimistic or pessimistic in those regions, leading the acquisition function astray.

Case study: 2D diffuser optimization with CFD

A diffuser is an expanding duct that reduces velocity and recovers pressure. We optimized one with 5 geometric parameters: inlet width, outlet width, length, expansion angle, and wall curvature. Two objectives (converted to single-objective via constraints): maximize pressure recovery coefficient Cp, subject to total pressure losses not exceeding a threshold. Each evaluation: OpenFOAM simpleFoam, mesh of approximately 200k cells, approximately 20 minutes on 8 cores.

Latin Hypercube: 15 initial points
GP with Matern 5/2 kernel
Acquisition: Expected Improvement
60 total evaluations (15 initial + 45 Bayesian iterations)

Result: optimum quality equivalent to NSGA-II at 500 evaluations. An 88% reduction in CFD computation time. This is the real power of Bayesian optimization: not that it finds better optima than genetic algorithms (it typically finds comparable optima), but that it finds them with an order of magnitude fewer function evaluations. When each evaluation is a CFD simulation, this translates directly into weeks of saved time.

The diffuser problem is particularly well-suited to Bayesian optimization because the response surface — pressure recovery as a function of geometry — is relatively smooth and has a clear global optimum. Problems with many local optima or discontinuities (e.g., flows with separation onset) are more challenging and may require larger initial designs or specialized kernels that encode knowledge about the physics.

Advantages over genetic algorithms

Sample efficiency: 10-100x fewer evaluations. Ideal when each evaluation is expensive. This is the single biggest advantage.
Uncertainty quantification: GPs naturally provide confidence intervals. You know not just what the predicted performance is, but how confident you should be in that prediction.
No population tuning: No need to adjust population size or genetic operators. Bayesian optimization has fewer hyperparameters and they are more robust to their settings.
Noise tolerant: Works with noisy observations, such as experimental data (not just simulation). The GP's noise parameter naturally handles measurement uncertainty.

Disadvantages: limited to approximately 20 design variables (GPs scale poorly with dimensionality — the number of points needed to maintain a given density grows exponentially), and sequential nature (harder to parallelize — each new evaluation depends on all previous ones). For high-dimensional problems, Bayesian optimization can be combined with dimensionality reduction techniques (PCA on design variables, active subspaces) or replaced with alternatives like random embeddings (REMBO) or trust-region Bayesian optimization (TuRBO).

Multi-objective Bayesian optimization

Several approaches extend Bayesian optimization to multiple objectives:

ParEGO: Scalarizes multiple objectives with random weights at each iteration, then optimizes with standard single-objective Bayesian optimization. The randomness ensures that over many iterations, the entire Pareto front is explored.
Expected Hypervolume Improvement (EHVI): Direct multi-objective acquisition function. Evaluates how much hypervolume a new point would add to the Pareto front. The gold standard but computationally expensive — requires integration over the objective space.
Multi-task GPs: Share information between related objectives through multi-task kernels, exploiting correlations between objectives. If two objectives are positively correlated, observations of one provide information about the other without additional evaluations.

In practice, ParEGO is the most accessible approach for problems with 2-4 objectives — it is simple to implement, computationally efficient, and the random weights ensure good coverage of the Pareto front. EHVI is preferable when the Pareto front has complex geometry (disconnected regions, sharp corners) that ParEGO's scalarization might miss. Multi-task GPs excel when objectives are known to be correlated — for example, maximizing lift and minimizing drag are usually anti-correlated, so learning one helps predict the other.

Practical recommendations

Matern 5/2 kernel as default for engineering (realistic smoothness for physical response surfaces). The RBF kernel is too smooth and leads to overconfident predictions.
Normalize inputs to [0, 1] and standardize outputs to zero mean and unit variance. This prevents numerical issues and makes kernel hyperparameter optimization more stable.
Multistart when optimizing the acquisition function: several random restarts avoid local optima of the surrogate itself. The acquisition function is often multi-modal even when the true objective is smooth.
Periodic validation of the GP: hold out a validation set to verify the surrogate is not drifting from reality. If the GP systematically overestimates or underestimates performance in certain regions, the acquisition function will be misled.
Nugget of 1e-6 on the kernel diagonal if the GP is poorly calibrated (underestimated uncertainty). Increase to 1e-3 or 1e-2 if numerical Cholesky failures occur.
Parallelization: Instead of a single recommendation, generate k candidates per iteration (q-EI, fantasizing) and launch several simulations in parallel. This maintains the sample efficiency advantage while exploiting available computing resources.

The most common failure mode in Bayesian optimization is an overconfident GP that predicts near-zero uncertainty in regions far from training data. This causes the acquisition function to ignore promising unexplored regions, converging prematurely to a suboptimal design. Diagnostics: plot the GP predictions across the design space (or through 1D slices for higher dimensions) and verify that uncertainty is large in regions without training points. If uncertainty collapses too quickly, increase the nugget or switch to a less smooth kernel like Matérn 5/2.

Bayesian optimization: surrogate models for expensive CFD