Cheap contouring of costly functions: the Pilot Approximation Trajectory algorithm

We consider a simple problem in which PAT is helpful: estimating the coefficients of a constant-coefficient steady state heat equation. The governing equation is

$$\begin{equation} a \nabla^2 u - bu = h, \end{equation} \tag{ 4 }$$

where a and b are the coefficients and h is a source term. This problem is two-dimensional in that the PDE depends on two coefficients; the domain of the PDE is three-dimensional. For example, the coefficients a and b might be parameters related to heat conductivity and perfusion; see Pennes (1948). We seek a 95% simultaneous confidence set for the unknown coefficients a and b based on observations y = (y₁,...,y_n)^T from u,

$$\begin{equation} y = G(u) + \epsilon, \end{equation} \tag{ 5 }$$

where is a n-vector representing observational noise and G = (G₁,...,G_n)^T represents the relationship between the observations and the underlying PDE. For example, G can be a linear interpolation operator which samples u of (4) in given measurement locations or it can be the negative temperature coefficient (NTC) resistor model. We assume that the components of are independent, identically distributed (iid) Gaussian variables with mean 0 and variance 1.

We will be contouring a function of $x=(a,b)\in \Delta \subset \mathbb R_+^2$, the vector of unknown parameters. Each value of x is a possible value for the coefficients in the heat equation, not a spatial location within the domain of the equation. Equation (5) has unique solution u for all positive x = (a,b) > 0. The solution can be approximated numerically, for example, using a finite-element method. Define the function f such that f(x) = G(u). That is, f maps a coefficient vector x = (a,b) into the n noise-free measurements according to the finite-element solution of (4) for coefficients x = (a,b). To evaluate f at a single point x = (a,b) requires finding the finite-element solution to equation (4) for the coefficients x = (a,b) and then applying the measurement model G to that finite-element solution.

The forward problem that maps a parameter vector to the observations thus can be written

$$\begin{equation} y = f(x) + \epsilon. \end{equation} \tag{ 6 }$$

We want to find a confidence set for x with confidence level 95%. This is a set C(y) for which

$$\begin{eqnarray*} &&\Pr \{ C(Y) \ni x \} \geqslant 95\% \end{eqnarray*}$$

whatever the true value x = (a,b) happens to be. An example of a 95% confidence set is

$$\begin{eqnarray*} &&C(y) = \left\{ x \in \Delta: \sum_{k=1}^n (y_k-f_k(x))^2 \leqslant \chi^2 \right\}, \end{eqnarray*}$$

where the constant χ is chosen such that

$$\begin{eqnarray*} &&\Pr \{ \| \epsilon \|^2 \leqslant \chi^2 \} \geqslant 0.95. \end{eqnarray*}$$

The appendix gives more explanation.

The boundary of the confidence set C(y) is the level-χ contour of the function g(x) =  $\left (\sum _{k=1}^n (y_k- f_k(x))^2\right )^{1/2}$. Evaluating g for a given parameter pair x = (a,b) can be expensive, especially if the problem is modeled in three spatial dimensions and the computational mesh contains large number of nodes. Regular contouring approaches might require evaluating g on a grid of 128 × 128 different parameter values, for instance. That would mean running the finite-element solver 128² times, once for each of the 128² parameter values. Clearly, it would be useful to be able to find the boundary of the confidence set with fewer evaluations of g; that is, the problem PAT addresses.

There are at least two ways to construct a pilot approximation $\hat g$ for PAT:

• (i)  
  Linear interpolation of g on a sparse grid. Evaluate g on a small grid of points (say 10 × 10 grid) in Δ. Divide each grid element into two triangles, forming a triangular tessellation of Δ. Set the pilot approximation $\hat g$ to be the piecewise linear function equal to g at the nodes, linear within each triangle. (Instead of linear interpolation, one might use MARS (Friedman 1991) or a Gaussian process model (O'Hagan and Kingman 1978, Sacks et al 1989, Kennedy and O'Hagan 2001)).
• (ii)  
  Solve the PDE at a low level of approximation for a dense grid of values of x. Computing g involves the finite-element solution of a PDE on a computational mesh. We can use the same finite-element solver on a coarser computational mesh to evaluate $\hat g$. For example, calculating g might involve a mesh with 128³ nodes, and calculating $\hat g$ might involve a mesh with 16³ nodes. In order to evaluate g(x), we use the mesh with 128³ nodes to solve (4) with coefficients x = (a,b), apply the observation sampling G to the solution, and calculate g(x), the sum of squared differences between the data and the sampled solution. To evaluate the approximation $\hat g(a,b)$, we process exactly same way except we numerically solve the PDE in the coarser mesh with 16³ nodes.

The first approach is more generic and can be used with all kind of problems. However, the second approach is more natural in problems arising from PDEs, and allows the approximation $\hat g$ to be constructed without evaluating g at all.

Suppose that g is a continuous function of two real variables. Let $\Delta \equiv [a_1,b_1]\times [a_2,b_2] \subset \mathbb R^2$. We seek the contour

$$\begin{eqnarray*} &&\gamma_{g, \tau} \equiv \{x \in \Delta: g(x)=\tau\}, \end{eqnarray*}$$

for a given level τ. An approximation of γ_g,τ is denoted $\hat {\gamma }_{g, \tau }$. Without loss of generality, we assume that τ = 0, since we can always contour g − τ instead of g. We suppress τ from the notation unless τ ≠ 0. The contour γ_g consists of a union of maximally connected pieces, indexed by i.

Standard algorithms for approximating γ_g start by dividing Δ into N smaller rectangles (cells) {Δ_k}^N_k=1. They also rely on knowing the value of g on the grid defined by the vertices of the cells. They approximate g at points not on the grid by interpolation, typically linear interpolation. Such algorithms rely on assumptions like the following:

• (i)  
  The function g is continuous.
• (ii)  
  The area of γ_g is zero.
• (iii)  
  Each piece of γ_g is either a closed curve or both of its endpoints are on the boundary of Δ.
• (iv)  
  The function g is equal to the level τ at most once on any edge of any cell.
• (v)  
  No piece of γ_g is entirely confined to the interior of a cell.

Assumptions (ii) and (iii) imply that g does not have a local extrema at which the value of g is 0. Assumptions (iv) and (v) ensure that the cells are small enough to find all the pieces of the contour. Aramini (1980) proposes approaches to obviate the need for assumption (iv) by adaptively subdividing cells; see also Suffern (1990). The method of Aramini (1980) requires knowing the approximate magnitude of the Laplacian of g, but we are avoiding the assumption that g is differentiable. Instead, we assume that the cells are sufficiently small that γ_g crosses no edge more than once. The PAT confidence bound on its own error—the area of the symmetric difference between the real and approximate lower level sets—does not depend on any of these assumptions.

Figure 1 sketches how a typical contouring algorithm works. First, the algorithm compares values of g at adjacent vertices until it finds an edge that the contour crosses. When it finds such an edge, the algorithm uses linear interpolation to estimate the point at which the contour intersects that edge. Then the algorithm moves to an adjacent cell, finds an edge of that cell crossed by the contour, and uses linear interpolation to estimate where the contour intersects this second edge. This continues until the algorithm either returns to the cell found in the first step, or reaches the boundary of the domain.

Zoom In Zoom Out Reset image size

If the algorithm returns to the cell found at the first step, this piece of the contour is a closed curve; the points where the contour was estimated to intersect the edges of cells are connected by straight lines (traced) to form a closed polygon that approximates the corresponding piece of γ_g. On the other hand, if the algorithm reaches the boundary of the domain Δ before returning to the first cell, assumption (iii) guarantees that this piece of the contour intersects the boundary of Δ twice. To find the other point where the contour intersects the boundary of Δ, the algorithm returns to the first point found and moves in the opposite direction until it reaches the boundary. Then the points at which this piece of the contour were estimated to intersect cell edges are joined by straight lines (traced) to form a piecewise linear curve that approximates the corresponding piece of γ_g.

The contour might consist of disjoint pieces, so when the algorithm finishes tracing one piece, it checks 'unvisited' cells to see whether the contour has another piece. If so, that piece is traced as described above. This process continues until every cell has been visited. If no edge is crossed by a contour, assumption (v) guarantees that the contour is the empty set.

As mentioned above, PAT is a modification of contouring algorithms like those just described. PAT assumes that g can be evaluated anywhere in Δ. PAT is designed for situations where evaluating g is expensive, so it seeks to keep the number of evaluations small but still approximate γ_g accurately.

PAT uses a pilot approximation $\hat g$ to determine where to evaluate g. For instance, if g involves an expensive simulation or physical experiment, $\hat g$ might be a spline function fitted to values of g at a small number of points in Δ. If computing g requires solving a PDE using a finite-element method, $\hat g$ might be based on a finite-element model using a coarse discretization.

In general, the contours of g and $\hat g$ intersect different cells. To find a cell intersected by γ_g, PAT starts at a point on $\gamma _{\hat g}$, the contour of $\hat g$. It searches a line perpendicular to $\gamma _{\hat g}$ for a cell intersected by the contour of g. (It should find such a cell eventually unless γ_g and $\gamma _{\hat g}$ are nearly orthogonal—which would imply that $\hat g$ is a poor approximation to g. In that case, PAT advances to the next edge intersected by $\gamma _{\hat g}$ and tries again.) Once such a cell is found, PAT uses the 'trajectory' (local orientation) of the contour of $\hat g$ to predict that of the contour of g. For instance, if the contour of $\hat g$ intersects the left and top edges of the cell in which the search started, PAT starts by assuming that the contour of g intersects the left and top edges of the cell just found by the search. If that assumption turns out to be false, PAT searches the other edges of the cell crossed by the contour of g to find an edge intersected by γ_g. If the assumption is true, PAT evaluates g fewer times, but if the assumption is false, PAT still works.

Like standard contouring algorithms, PAT divides the domain into cells {Δ_k}^N_k=1. PAT requires both g and $\hat g$ to satisfy assumptions (i)–(v). To simplify the exposition, we assume in this section that the contours of g and $\hat g$ have same number of disjoint pieces; we relax this assumption in section 2.4. And we assume in this section that γ_g and $\gamma _{\hat g}$ each intersect at most two edges of any cell; we relax this assumption in section 2.3.

Below are the steps in the basic PAT algorithm. The cells intersected by the ith piece of the approximate contour are denoted $\{\hat \Delta ^{i}_{j'}\}$, j' = 0,...,J'_i. Cells intersected by the ith piece of γ_g are denoted {Δⁱ_j}, j = 0,1,...,J_i.

• 0.  
  Construct the pilot approximation $\hat g$ of g. Divide the domain Δ into N cells, {Δ_k}^N_k=1.⁴ Set k = i = 1.
• 1.  
  Mark cell Δ_k as visited.
• 2.  
  • (a)  
    If k = N and Δ_k has no edge for which $\hat g(x_1)\hat g(x_2)\leqslant 0$, where x₁ and x₂ are the vertices of the edge, exit.
  • (b)  
    Otherwise if k < N and Δ_k has no edge for which $\hat g(x_1)\hat g(x_2)\leqslant 0$ , increment k and return to step 1.
  • (c)  
    Otherwise Δ_k has an edge with $\hat g(x_1)\hat g(x_2)\leqslant 0$, so $\gamma _{\hat g}$ intersects Δ_k. The assumptions guarantee that it has at least two such edges. Denote two of them $\hat e^{i}_{0}$ and $\hat e^{i}_{-1}$. Set $\hat \Delta ^{i}_{0}=\Delta _k$.
• 3.  
  Find an edge intersected by γ_g as follows:

  • (a)  
    Use linear interpolation to estimate where $\gamma _{\hat g}$ intersects the edges $\hat e^{i}_{0}$ and $\hat e^{i}_{-1}$. Denote those points by $\hat {x}_{0}^{i}$ and $\hat {x}_{-1}^{i}$. If $\hat g$ has a zero at a vertex of the cell, we can choose that vertex to be one of the points.
  • (b)  
    Find the line perpendicular to the line segment between the points $\hat {x}_{0}^{i}$ and $\hat {x}_{-1}^{i}$ that goes through $(\hat {x}_{0}^{i}+\hat {x}_{-1}^{i})/2$; see figure 2. Search for a root of g on this line that is not on a piece of the contour already traced (our implementation of PAT uses Müller's method (Müller 1956); other methods might be more efficient). If no new root is found on this line, increment k and return to step 1. If a root is found, denote it xⁱ₀. The point xⁱ₀ is on the piece of γ_g we will trace next.
  • (c)  
    Find the cell containing xⁱ₀ and denote it Δⁱ₀.
• 4.  
  Follow and trace the piece of γ_g just found: set j = j' = 0. Then:

  • (a)  
    Predict which edge γ_g crosses next. The prediction is denoted e_c. If $\gamma _{\hat g}$ is in the same cell (i.e., has also entered the cell Δⁱ_j+1 through the edge eⁱ_j or, on the first step, started from the same cell), set e_c to be an edge that $\gamma _{\hat g}$ crosses next. Otherwise, let e_c be the edge next crossed by the function $\hat g(\cdot )-\hat g(x_{j}^{i})$, the pilot approximation $\hat g$ shifted by its value at the last contour point of g, xⁱ_j.
  • (b)  
    If edges e_c and eⁱ_j are adjacent, they share a vertex at which g is already known. In that case, evaluate g at the other vertex of e_c to check whether γ_g crosses e_c. If it does, set eⁱ_j+1 = e_c.
  • (c)  
    If e_c is parallel to eⁱ_j, evaluate g at one vertex of e_c (e.g. the vertex with the smaller absolute value of $\hat g$). Then g is known at both vertices of one edge perpendicular to eⁱ_j; check whether γ_g crosses that edge. If so, that edge is eⁱ_j+1. If not, evaluate g at the remaining vertex of Δⁱ_j+1 to find the edge eⁱ_j+1 that γ_g intersects.
  • (d)  
    Find the root xⁱ_j+1 of g on edge eⁱ_j+1.
  • (e)  
    If $\gamma _{\hat g}$ also crosses edge eⁱ_j+1, trace the portion of $\gamma _{\hat g}$ between the edge $\hat {e}^{i}_{j'}$ and the cell Δⁱ_j+1 the same way the standard algorithm does. Update the index j' and mark the cells $\hat {\Delta }^{i}_{j'}$ intersected by that portion of $\gamma _{\hat g}$ as visited.
  • (f)  
    Move to the cell that is across edge eⁱ_j+1 and denote it by Δⁱ_j+1.Set j←j + 1 and repeat these steps until either the initial cell Δⁱ₀ or the boundary is reached. If a boundary is reached, return to Δⁱ₀ and trace in the other direction, as the standard contour algorithm does. If the end of $\gamma _{\hat g}$ has not been reached, trace $\gamma _{\hat g}$ to its end and mark all cells on its path as visited. Increment i.
• 5.  
  If there are cells not yet visited, increment k until Δ_k is not a visited cell and continue from step 1.
• 6.  
  For each i, connect the points {xⁱ₀,xⁱ₁,...,xⁱ_Ji} with straight line segments to form the approximation $\hat {\gamma }_{g}$ of γ_g.

Zoom In Zoom Out Reset image size

To approximate contours at additional levels, the pilot approximation $\hat g$ could be reused, or $\hat g$ could be improved by incorporating all the points at which g has been evaluated.

In the previous section we assumed that the contour intersects at most two edges of any cell. This may not hold if (a) the cell contains a saddle point (the contour crosses itself or two pieces of the contour intersect the same cell) or (b) there is a local minimum or maximum in the cell (so that the same piece of the contour intersects the cell four times, or at three edges and a vertex).

There have been attempts to distinguish these situations automatically, for example by using linear interpolation to decide which edges to connect (Aramini 1980) or by evaluating the function at the center of the cell and comparing that value to the values at vertices (Aramini 1980). Instead, we follow the contour in small steps from one edge to another edge. We believe that this approach errs less frequently in practice, but we have not characterized the conditions under which that is true.

Suppose PAT is at step 4(a) and it turns out that γ_g crosses all four edges of the cell. Then the following approach is used to determine the next edge eⁱ_j+1 and the point xⁱ_j+1. The idea is to follow the contour from xⁱ_j to another edge. One way to do this is to move a short distance along the tangent to γ_g and then move perpendicular to that tangent to get back to γ_g, and repeat these steps until we cross an edge. However, since g is not necessarily differentiable, the tangent to γ_g might not exist. Hence, we modify this prediction–correction approach.

We construct a step direction as follows. Pick a small value δ > 0 (for instance, 10% of length of edge eⁱ_j). Find the unit normal n to edge eⁱ_j. Set k←0, d₀←δn, and y₀←xⁱ_j. While y_k∈Δⁱ_j+1:

• 1.  
  Set $\hat y_{k+1} \leftarrow y_k +d_k$.
• 2.  
  On the line through $\hat y_{k+1}$ perpendicular to d_k, find the closest root of g. Set y_k+1 equal to that root.
• 3.  
  Set d_k+1 = δ(y_k+1 − y_k) and let k←k + 1.

This approach is illustrated in figure 3.

Zoom In Zoom Out Reset image size

The edge eⁱ_j+1 is chosen to be the edge crossed by the line segment between y_k and y_k+1. The point xⁱ_j+1 is chosen to be the intersection point on this line segment. The points y₁,...,y_k can be included in the series of points connected to form the approximation $\hat {\gamma }_{g}$ to the ith piece γ_g.

So far, we have assumed that γ_g and $\gamma _{\hat g}$ have the same number of pieces. If $\gamma _{\hat g}$ has more pieces than γ_g, PAT will not succeed in finding a new piece of γ_g when tracing the extra pieces of γ_g, potentially wasting a large number of evaluations of g—but otherwise harmless. On the other hand, if γ_g has more pieces than $\gamma _{\hat g}$, PAT might not find them all. This can matter in applications. For instance, missing a piece of γ_g can cause a misfit-based confidence set to have a confidence level less than its nominal confidence level, by erroneously excluding parameter values that belong in the set. This section discusses how to augment PAT to make it more likely to find all the pieces of γ_g and how to measure statistically the area of the symmetric difference between the lower level set of g and the estimated lower level set of g, as discussed in section 1.

2.4.1. Heuristic approaches.

2.4.2. A priori knowledge about g.

2.4.3. A priori knowledge about $\vert \hat g - g \vert$.

If we have a bound on the error of $\hat g$ as an approximation to g, we can use it to check whether PAT has missed pieces of γ_g. For example, if we knew a function M(x) for which

$$\begin{eqnarray*} &&\vert g(x)-\hat g(x)\vert \leqslant M(x)\quad\textrm{for all } x\in\mathcal{X}, \end{eqnarray*}$$

it would suffice to evaluate g at grid points x_j for which $\vert \hat g(x_j)\vert \leqslant M(x_j)$ (perhaps starting with points for which $\vert \hat g(x_j)\vert M(x_j)^{-1}$ is smallest) and check the sign of g at those points.

2.4.4. Random sampling and contour error.

We can use random sampling both to look for missing pieces and to collect statistical evidence about the accuracy of the PAT contour $\hat {\gamma }_g$ as an approximation to γ_g. Recall from section 1 that the difference between the true lower level set of g at level τ, Λ_τ, and the estimated lower level set of g, $\hat {\Lambda }_\tau$, is the set on which the approximate contour $\hat {\gamma }_g$ errs in classifying points in Δ. As discussed in the introduction, the approximate contour $\hat {\gamma }_g$ makes classification errors on two sets: points where the true value of g is strictly positive but the approximate contour claims that g is nonpositive, and points where the true value of g is nonpositive but the approximate contour claims that g is positive. The area of the union of those sets—the symmetric difference $\Lambda _\tau \ominus \hat {\Lambda }_\tau$ of Λ_τ and $\hat {\Lambda }_\tau$—is a measure of the accuracy of the approximate contour. Normalizing that area by the total area of Δ gives the fraction of the area of the domain Δ misclassified by the approximate contour.

To keep things simple, we consider how to estimate the fraction p of the area of Δ that is misclassified by the PAT contour $\hat {\gamma }_g$—the relative area of the symmetric difference. The chance that a point selected at random uniformly from Δ is misclassified by $\hat {\gamma }_g$ is p. If a point near $\hat {\gamma }_g$ is misclassified, it could mean that $\hat {\gamma }_g$ is slightly out of place. If a point far from $\hat {\gamma }_g$ is misclassified, that suggests that $\hat {\gamma }_g$ is a missing piece of γ_g. To check, we can search for a root of g on any line segment between that point and the nearest piece of $\hat {\gamma }_g$. If we find such a root, it becomes a starting point to trace a piece of γ_g that had been missed.

On the other hand, suppose we draw n points at random, independently, uniformly from Δ, and the sign of g at every point in the sample agrees with its sign according to the PAT contour $\hat {\gamma }_g$: no point in the sample is in the symmetric difference $\hat {\Lambda }_\tau \ominus \Lambda _\tau$. Then an event has occurred that has the probability (1 − p)ⁿ. We can use this to find a 1 − α upper confidence bound on p: for this event to have a probability greater than α, p must be less than 1 − α^1/n. For instance, for n = 50, we would have 95% confidence that P < 5.8%; for n = 100, we would have 95% confidence that P < 3%.

More generally, if n independent draws find m misclassified points, a 1 − α upper confidence bound for the relative area of $\hat {\Lambda }_\tau \ominus \Lambda _\tau$, the part of the domain Δ misclassified by the PAT contour, is

$$\begin{equation} u(n, m, \alpha) \equiv \max \left \{ p: \sum_{i=0}^m p^i (1-p)^{n-i} \geqslant \alpha \right \}. \end{equation} \tag{ 7 }$$

Evaluating g at a random sample of points in Δ allows us to have confidence in the accuracy of the PAT contour, despite the fact that the assumptions behind PAT might be wrong—at the expense of n extra evaluations of g.

In some applications, one kind of error is more important than the other. For instance, erroneously including extra points in a confidence set makes the set conservative, while erroneously excluding points from a confidence set is misleadingly optimistic. Moreover, the 'cost' of misclassifying points might be higher for some points than for others. For instance, erroneously excluding from a confidence set of climate models a model in which global temperature rises 10 K in 20 years might be a more serious error than erroneously excluding a model in which global temperature rises 0.5 K in 20 years. Such differences can be taken into account by using weighted areas (integration against a nonuniform measure instead of Lebesgue measure) and carrying out the above sampling using corresponding nonuniform distributions.

Let g be a sum of 2 two-dimensional Gaussian bumps:

$$\begin{eqnarray*} &&g(x)=\mathrm{e}^{-(Rx-\mu_-)^{\mathrm{T}} C (Rx-\mu_-)}+\mathrm{e}^{-(Rx-\mu_+)^{\mathrm{T}} C (Rx-\mu_+)},\quad x\in\mathbb R^2 , \end{eqnarray*}$$

where the matrices R and C and the vectors μ_± are given by

$$\begin{eqnarray*} R=\left(\begin{array}{@{}cc@{}} \cos(\frac\pi 4) & \sin(\frac\pi 4) \\[6pt] -\sin(\frac\pi 4) & \cos(\frac\pi 4) \end{array} \right),\quad C=\left(\begin{array}{@{}cc@{}} 8 & 0 \\ 0 & 10 \end{array} \right), \quad \mu_\pm=\left(\begin{array}{@{}c@{}} \pm0.6\\0\end{array}\right) . \end{eqnarray*}$$

The matrix R rotates the coordinate system 45° clockwise. The domain for the contour plot is Δ ≡ [ − 1,1]². The function g is plotted in figure 4.

Zoom In Zoom Out Reset image size

This function has a saddle point at the origin. To test the ability of PAT to contour saddle points, we contour at level τ = g(0,0), the value of g at the saddle point.

The pilot approximation $\hat g$ was chosen to be the piecewise linear approximation of g using an equispaced 9 × 9 grid in Δ (see description in section 1) and it is plotted in figure 4.

The PAT contour $\hat {\gamma }_g$ found using a 32 × 32 grid is plotted in figure 5. For comparison, figure 5 shows a contour computed using the built-in contouring algorithm in R with data on a 1024 × 1024 grid and on the 32 × 32 grid. The 1024 × 1024 grid is dense enough that the approximate contour is essentially the true contour of g. Figure 6 is a close-up of the contour near the saddle point.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The PAT contour is almost indistinguishable from the standard contour on a 1024 × 1024 grid. It is comparable to the standard contour on a 32 × 32 grid except near the saddle point, where the PAT contour is more accurate. PAT correctly finds that the contour crosses itself, even though the contour of the pilot approximation does not. Moreover, PAT is more accurate than the standard contouring algorithm on the 32 × 32 grid, in part because PAT locates zeros of g with a root finder rather than approximating them by linear interpolation. Indeed, if we treat the contour on the 1024 × 1024 grid as ground truth, the standard contouring on a 32 × 32 grid misclassifies 0.38% of the domain, while PAT misclassifies only 0.12% of the domain (the area of the symmetric difference was computed using the gpclib library in R. Repeated random samples of 100 points generally found none that were misclassified by the PAT contour (no sample point was in $\hat {\Lambda }_\tau \ominus \Lambda _\tau$). The 95% upper confidence bound on the fraction of Δ misclassified by the PAT contour based on a sample of size 100 that shows no misclassifications is about 3%, rather larger than the true percentage, 0.12%. Occasionally, a random sample of 100 points revealed one point misclassified by the PAT contour; the corresponding upper 95% confidence bound is 4.7%. By sampling more than 100 points, these upper confidence bounds could be reduced. However, the bounds decrease roughly like the reciprocal of the square root of the sample size. For instance, to reduce the 95% upper confidence bound from 3% to 0.2% would require evaluating g at least 1500 times, to reduce the sampling variability in the accuracy assessment—not to improve the contour.

To compute the PAT contour required evaluating g 370 times (the pilot approximation $\hat g$ was evaluated 1109 times). To build the pilot approximation $\hat g$, g was evaluated 9² = 81 times, so g was evaluated 451 times in all. For comparison, computing the contour on a 32 × 32 grid using a standard contouring algorithm requires 32² = 1024 evaluations of g. In this example, our algorithm requires 56% fewer evaluations of g. Including the 100 evaluations of g at randomly selected points, the total number of evaluations is still 46% less than the standard contouring algorithm requires—and gives a more accurate contour and an internal measure of accuracy, which the standard method lacks. (Of course, the same sampling approach could be used to assess the accuracy of the standard method, at the cost of an additional 100 evaluations of g for parity.) The number of evaluations of g that PAT requires depends on the root-finding algorithm and its tuning constants; we used the algorithm of Müller (1956).

Missing contour pieces.

To test the ability of PAT to find a contour accurately even when the topology of the contour of the pilot approximation is wrong, we shifted the level τ for which the contour is constructed by 0.005 (i.e. τ = g(0,0) + 0.005). This makes the contour γ_g,τ consist of two closed curves. We also shifted the approximation $\hat g$ by 0.045 so that the contour of the pilot approximation $\gamma _{\hat g, \tau }$ is a single closed curve that encircles both pieces of the contour γ_g,τ, so the contour γ_g,τ has more pieces than $\gamma _{\hat g, \tau }$, the contour of the pilot approximation.

The method described in section 2.4 finds the piece. That method checks for whole-cell 'gaps' between the PAT contour $\hat {\gamma }_{g, \tau }$ and the approximate contour of the pilot approximation, $\hat {\gamma }_{\hat g, \tau }$. In this example, the missing piece of the PAT contour is inside the approximate contour, so the gaps are outside the PAT contour. Since g is below τ near but outside the PAT contour, the algorithm will detect that a piece of the contour is missing if it finds a point x outside the traced piece for which g(x) > τ. To look for such a point, it evaluates g within the gap at the grid point where $\hat g$ is largest, since, if $\hat g$ is a good approximation of g, we expect g to be large there, too. If g > τ at that point, the algorithm searches the line segment between that point and the closest point at which g has been already evaluated (a vertex of the cell intersected by the contour) to find a point where g = τ. It then traces the (new) piece of the contour that passes through that root of g − τ.

The PAT contour on a 32 × 32 grid is shown in figure 7. Computing the PAT contour required 423 evaluations of g: 8² = 64 to find the pilot approximation and 359 to trace the contour. If we consider the contour on the 1024 × 1024 grid to be the ground truth, the standard contour on a 32 × 32 grid misclassifies 0.33% of the domain, while PAT misclassifies only 0.15% of the domain. Again, PAT gives a more accurate contour than the standard method, using far fewer evaluations of g.

Zoom In Zoom Out Reset image size

Workload scaling.

The steady-state temperature distribution T(·) of the gas is modeled by the heat equation. Since there are no internal sources, this equation is:

$$\begin{equation} \rho C\left(v\cdot\nabla T\right)-\nabla\cdot (\kappa\nabla T)= 0, \end{equation} \tag{ 8 }$$

where C is heat capacity, v is the velocity field and κ is thermal conductivity. The velocity field is modeled by the Navier–Stokes equation:

$$\begin{equation} \rho(v\cdot\nabla)v-\nabla\cdot(2\mu\epsilon)+\nabla p=\rho h, \end{equation} \tag{ 9 }$$

where ρ is density, μ is viscosity, p is pressure, h is gravitational force and is the infinitesimal strain tensor. The velocity and temperature boundary conditions at the input boundary vary quadratically with height. The velocities at the walls are zero and purely horizontal at the output boundary. The temperature at the walls is fixed, and the flow satisfies the Neumann zero boundary condition at Γ_out.

We solved the heat equation and the Navier–Stokes equation as coupled PDEs using Elmer. The function f and its pilot approximation $\skew5\hat{f}$ (see the notation in the appendix) are as follows. The function f maps (x₁,x₂) to HT, where T is the temperature distribution for a fixed viscosity of 10^−x₁ Ns m⁻² and a fixed thermal conductivity of 10⁻²·x₂ W (m K)⁻¹, and H is the matrix that interpolates linearly between points on the computational mesh to calculate T at the measurement points. The domain of f is Δ = [3,4] × [2,3], which we verified experimentally to contain every point in the confidence set. A logarithmic transformation of viscosity improves numerical stability, since the effect of viscosity on velocity and temperature is exponential. We scaled the x₂-axis to make the axes equally long to improve computational efficiency. An approximation $\skew5\hat{f}$ for f is formed in the same way using a very sparse mesh.

We used PAT to approximate the contour corresponding to the boundary of the 95% χ² confidence set (see the appendix). The function g was defined using (A.2) with p = 2 (χ² confidence sets). The pilot approximation $\hat g$ is obtained similarly by replacing f with $\skew5\hat{f}$ in equation (A.2).

The computational mesh used to evaluate f contained 330 elements and 1099 nodes. For the pilot approximation $\skew5\hat{f}$, the mesh contained 108 elements and 389 nodes. Elmer becomes unstable in this problem if there are significantly fewer elements than that. The material parameters used in the computations are shown in table 2. PAT contours were computed using a grid of 30 × 30 cells (31 × 31 grid points).

The simulated data were computed using Elmer with a dense mesh of 563 elements and 1830 nodes. Independent, identically distributed zero-mean Gaussian noise with standard deviation 2.79 K (2% of the elevation in the temperature) was added to the simulated values.

To assess accuracy of the PAT confidence sets, we also computed the contour using the standard algorithm in R on a 256 × 256 grid. The PAT contours and reference contours are shown in figure 9. The PAT contours are almost indistinguishable from contours computed by standard methods on the dense grid. Differences are visible only in sharp corners where the cell size is too large for a straight line segment to turn as rapidly as the actual contour.

Zoom In Zoom Out Reset image size

Computing the PAT contour for the χ² confidence set required evaluating the function f 243 times and the function $\skew5\hat{f}$ 1023 times. The number of the evaluations of f was about a quarter of 31² = 961, the minimum number of evaluations a standard contouring method would use to compute the contour on a 31 × 31 grid. Moreover, the PAT contour is more accurate than the standard contouring method on a 31 × 31 grid: the relative area of $\hat {\Lambda }_\tau \ominus \Lambda _\tau$ for the PAT confidence set is 0.06% while the relative area for a standard contour on the same grid is 0.15%.