The Algorithm of Gu and Eisenstat and D-Optimal Design of Experiments

Abstract

This paper addresses the following problem: given m potential observations to determine n parameters, $m>n$, what is the best choice of n observations. The problem can be formulated as finding the $n\times n$ submatrix of the complete $m\times n$ observation matrix that has maximum determinant. An algorithm by Gu and Eisenstat for a determining a strongly rank-revealing QR factorisation of a matrix can be adapted to address this latter formulation. The algorithm starts with an initial selection of n rows of the observation matrix and then performs a sequence of row interchanges, with the determinant of the current submatrix strictly increasing at each step until no further improvement can be made. The algorithm implements rank-one updating strategies, which leads to a compact and efficient algorithm. The algorithm does not necessarily determine the global optimum but provides a practical approach to designing an effective measurement strategy. In this paper, we describe how the Gu–Eisenstat algorithm can be adapted to address the problem of optimal experimental design and used with the QR algorithm with column pivoting to provide effective designs. We also describe implementations of sequential algorithms to add further measurements that optimise the information gain at each step. We illustrate performance on several metrology examples.

1. Introduction

While machine learning paradigms are based on access to large sets of data to extract the information of interest, many practical problems in science and technology arise in a context in which data are sparse and potentially expensive to generate. Given limited experimental resources, it is important that these resources are used effectively. The design of experiment (DoE) is an important area of study in applied mathematics and statistics. Many DoE studies relate to problems in which a response variable of interest, e.g., the yield of a crop, is possibly influenced by a number of control or stimulus variables, e.g., crop variety, soil type, irrigation regime, fertiliser treatment, etc. The form of the response function is usually unknown and the DoE problem is to design a number of experiments with the control variables set to different values, in order to estimate the response function from a set of candidate response functions and consequently choose control values that optimise the response; see, e.g., [1,2,3,4,5,6,7,8]. The DoE problem has been addressed in applications across many scientific disciplines, including the placement of sensors in a electrical power grid [9,10] and the measurement of thermal properties of materials [11,12], for example.

In this paper, we look at a particular aspect of the design problem: given a known response model, how do we choose a minimal number of measurements, taking into account the associated uncertainties, in order to estimate all the parameters of the system under study as accurately as possible. Thus, the problem can be posed as: given m potential observations to determine n parameters, $m>n$, what is the best choice of n observations. In the context of (linear) least squares estimation, a version of the problem can be formulated as finding the $n\times n$ submatrix of the complete $m\times n$ observation matrix that has maximum determinant. The problem can be posed as an example of the following general optimisation problem. Suppose $A=\{a_i,i=1,\dots ,m\}$ is a set with m elements. Represent the set $P\left(A\right)$ of all subsets of A as the set of characteristic functions $\chi :A⟶\{0,1\}$: if $B\subset A$, then the corresponding function ${\chi }_B\in P\left(A\right)$ is such that ${\chi }_B\left(a_i\right)=1$ if $a_i\in B$ and is 0 otherwise. Let $D:P\left(A\right)⟶\mathbb{R}$ be a real valued function defined on $P\left(A\right)$. The optimisation problem is

$$\underset{\chi \in P\left(A\right)}{min}D\left(\chi \right)\mathrm{subject}\mathrm{to}\sum _{i=1}^m\chi \left(a_i\right)=n.$$

(1)

In theory, the global solution to this DoE problem can be found by a brute force assessment of all choices of subsets, a total of $m!/(n!(m-n)!)$ possibilities. For all but small m and n, this approach will be computationally infeasible, even if the functional D is cheap to evaluate. Two other related approaches are commonly used to partially solve this type of combinatorial optimisation problem.

The first approach, with a long history, see e.g., [13,14,15], is to use a sequential approach that at each stage constructs a new candidate subset by removing one element from the subset, replacing it by a new element that leads to a strict improvement in D. The computational cost per step for such an algorithm will be given by the cost evaluating D and the cost of searching for a suitable exchange from the $n(m-n)$ possibilities. Ideally, the computational cost of an exchange should be no worse than $O\left(m\right)$ or $O(mlogm)$ as a function of m. The second approach is based on convex relaxation, replacing a combinatorial problem by a convex optimisation problem [16]. Let $0\le w_i=w\left(a_i\right)\le 1$, where we think of $w_i$ as specifying the degree to which the ith element is in the solution optimal subset. If $W\left(A\right)$ is the set of functions $A⟶[0,1]$, then the optimisation problem is defined by

$$\underset{w\in W\left(A\right)}{min}D\left(w\right)\mathrm{subject}\mathrm{to}\sum _{i=1}^{m}w\left(a_i\right)=n.$$

(2)

In this approach, the combinatorial constraint $\chi \left(a_i\right)\in \{0,1\}$ is relaxed to the convex constraint $0\le w\left(a_i\right)\le 1$. Often, the optimal w is in a fact a characteristic function associated with an n-element subset. The optimisation algorithms for such problems are generally iterative and often also have a sequential element, such as in sequential quadratic programming [17] and semidefinite programming [18]. Convex relaxation approaches to DoE problems were considered in [4,9,10,16], for example.

In this paper, we show how two algorithms, both following a sequential approach and developed for numerical linear algebra applications, can be repurposed to address the design problem of choosing n observations from m possible observations.

The remainder of this paper is organised as follows. In Section 2, we give an overview of linear least squares estimation, while in Section 3, we describe commonly used aggregate measures of uncertainty and the corresponding criteria used in optimal experimental design. In Section 4, we describe two algorithms that address the DoE problem under consideration, as well as sequential algorithms for choosing additional observations in a step-wise optimal way. Example applications of the algorithms to design problems arising in metrology  [19,20,21,22] are given in Section 5. Our concluding remarks are given in Section 6.

2. Least Squares Problems

For $m\times n$ observation matrix C, $m\times 1$ vector of measurements y, the least squares estimate a is the solution of

$$\underset{\mathsf{\alpha }}{min}{(\mathit{y}-C\mathsf{\alpha })}^{\top }(\mathit{y}-C\mathsf{\alpha }).$$

The QR factorisation of C expresses $C=QR$ as a product of an $m\times m$ orthogonal matrix Q and an upper triangular matrix R. If

$$C=\left[Q_1{Q}_2\right]\left[\begin{array}{c}{R}_1\\ \mathbf{0}\end{array}\right]=Q_1{R}_1,$$

where $Q_1$ is the submatrix constructed from the first n columns of Q and $Q_2$ that is constructed from the remaining $m-n$ columns, then

$$\mathit{a}={\left(C^{\top }C\right)}^{-1}{C}^{\top }\mathit{y}=R_1^{-1}{Q}_1^{\top }\mathit{y},$$

$$\widehat{\mathit{y}}=C\mathit{a}=C{\left(C^{\top }C\right)}^{-1}{C}^{\top }\mathit{y}=Q_1{Q}_1^{\top }\mathit{y},$$

$$r=\mathit{y}-C\mathit{a}=(I-Q_1{Q}_1^{\top })\mathit{y}=Q_2{Q}_2^{\top }\mathit{y},$$

are the least squares solution, model approximant, and associated residual vector. The QR factorisation enables a to be given by the solution of $R_1\mathit{a}=Q_1^{\top }\mathit{y}$ and this solution can be constructed using back substitution, solving for $a_n$ first, then $a_{n-1}$, etc., exploiting the upper-triangular form of $R_1$. The QR factorisation also avoids the numerical loss of accuracy that can occur in forming $C^{\top }C$ and its inverse in the case in which C is ill-conditioned [23] (Chapter 5). The computational cost of evaluation the least-squares solution in $O\left(m{n}^2\right)$.

The $n\times m$ matrix $S=R_1^{-1}{Q}_1^{\top }$ is the sensitivity matrix of the solution parameters a with respect to the data y with $s_{ji}=\partial a_j/\partial y_i$. If the variance matrix associated with y is $V_{\mathit{y}}$, then the variance matrix associated with the solution a is

$$V_{\mathit{a}}=S{V}_{\mathit{y}}{S}^{\top }.$$

For the case $V_{\mathit{y}}=I$,

$$V_{\mathit{a}}=S{S}^{\top }={\left(R_1^{\top }{R}_1\right)}^{-1}={\left(C^{\top }C\right)}^{-1}.$$

For this case, the trace of the variance matrix is given by

$$\mathrm{Tr}\left(V_{\mathit{a}}\right)=\sum _j\sum _i{s}_{ji}^2,$$

the variance associated with $a_j$ is $u^2\left(a_j\right)={\sum }_i{s}_{ji}^2$ and ${\sum }_j{s}_{ji}^2$ is the contribution of the uncertainty associated with the ith observation to the trace of $V_{\mathit{a}}$.

The least-squares solution a is the maximum likelihood estimate of $\mathsf{\alpha }$ for the model $\mathit{y}\in \mathcal{N}(C\mathsf{\alpha },I)$. If y is associated with the variance matrix $V_{\mathit{y}}$, i.e., $\mathit{y}\in \mathcal{N}(C\mathsf{\alpha },V_{\mathit{y}})$, and $V_{\mathit{y}}$ is full rank and has Cholesky decomposition $V_{\mathit{y}}=L{L}^{\top }$, where L is lower triangular [23], then the maximum likelihood estimate a of $\mathsf{\alpha }$ solves the Gauss–Markov problem

$$\underset{\mathsf{\alpha }}{min}{(\mathit{y}-C\mathsf{\alpha })}^{\top }{V}_{\mathit{y}}^{-1}(\mathit{y}-C\mathsf{\alpha }).$$

(3)

If

$$\tilde{C}=L^{-1}C,\tilde{\mathit{y}}=L^{-1}\mathit{y},$$

then the solution of (3) is the least-squares solution of $\tilde{C}\mathsf{\alpha }\approx \tilde{\mathit{y}}$:

$$\mathit{a}={\left({\tilde{C}}^{\top }\tilde{C}\right)}^{-1}{\tilde{C}}^{\top }\tilde{\mathit{y}}.$$

In particular, if $V_{\mathit{y}}$ is a diagonal matrix with $u_i^2>0$ in the ith diagonal element, then $\tilde{C}$ is constructed from the weighted rows of C, and similarly for $\tilde{\mathit{y}}$:

$${\tilde{\mathit{c}}}_i=w_i{\mathit{c}}_i,{\tilde{y}}_i=w_i{y}_i,w_i=1/u_i.$$

In general, the cost of evaluating the Cholesky factorisation of $V_{\mathit{y}}$ is $O\left(m^3\right)$, which could be problematic for large m. Often, $V_{\mathit{y}}$ has a structure that enables the factorisation to be performed in $O\left(m\right)$, as in the simple but common case of a diagonal variance matrix.

3. Aggregate Measures of Uncertainty

A central issue in design of experiment analysis is to maximise some measure of the information gain subject to constraints on resources. In a metrology context, in particular, it is natural to express the information gain in terms of uncertainties associated with the fitted parameters; i.e., related to the $n\times n$ variance matrix $V_{\mathit{a}}$. Given an aggregate measure of uncertainty derived from $V_{\mathit{a}}$, one experiment will be judged to be more informative than another if it is associated with a lower aggregate uncertainty. There are a number of aggregate measures that are commonly used. In this paper, we are concerned with two: A-optimality and D-optimality.

A full rank variance matrix V is a symmetric positive definite matrix and has an eigenvalue decomposition with eigenvalues ${\lambda }_j$, $j=1,\dots n$. The A-measure is the trace of $\mathrm{Tr}\left(V\right)$ of V, the sum of the diagonal elements of V, as well as the sum of the eigenvalues of V:

$$\mathrm{Tr}\left(V\right)=\sum _{j=1}^n{\lambda }_j.$$

We note that $\mathrm{Tr}\left({\sigma }^{2}V\right)={\sigma }^2\mathrm{Tr}\left(V\right)$. The D-measure is the determinant $|V|$ of V, and also the product of the eigenvalues of V:

$$|V|=\prod _{j=1}^n{\lambda }_j.$$

We also evaluate

$$\overline{d}\left(V\right)={|V|}^{1/n}={\left(\prod _{j=1}^n{\lambda }_j\right)}^{1/n},$$

(4)

the geometric mean of the eigenvalues, noting that $\overline{d}\left({\sigma }^{2}V\right)={\sigma }^2\overline{d}\left(V\right)$.

Arguments can be made for and against using each of these two measures. Given V, the A-measure can be evaluated in $O\left(n\right)$ steps, while the D-measure will usually involve $O\left(n^3\right)$ steps, using the Cholesky factorisation of the variance matrix to evaluate its determinant. The A-measure has the apparent advantage of being most easily interpreted as an aggregate uncertainty, being the sum of the variances ${\sum }_j{u}^2\left(a_j\right)$ associated with each of the parameter estimates $a_j$ or, more precisely, the numerical values of these estimates. From a dimensional analysis point of view [24], the A-measure can only be meaningful if all the parameters in a have the same dimension, length, mass, etc. Otherwise, the trace involves adding quantities of different dimensions. The D-measure is consistent with a dimensional analysis.

The A-measure can also be regarded as the variance associated with the sum ${\sum }_j{a}_j$ of parameter estimates, assuming these estimates as statistically independent. In general, they will not be independent, and one disadvantage of the A-measure is that it takes no account of correlation, being defined solely in terms of the diagonal elements of V. The D-measure does take into account correlation. Consider the $2\times 2$ variance matrix

$$V=\left[\begin{array}{cc}1& \rho \\ \rho & 1\end{array}\right],\mathrm{Tr}\left(V\right)=2,|V|=1-{\rho }^2,$$

associated with parameters $a_1$ and $a_2$. The A-measure is independent of the correlation coefficient $\rho$ while the D-measure depends directly on ${\rho }^2$. If the result of an experiment is to change $\rho$ from 0 to near 1, there is a gain of information, since $a_1-a_2$ is now much more accurately estimated. The D-measure reflects this information gain, the A-measure does not.

The D-measure is invariant to re-parametrizations of the model, while the A-measure is not. Suppose $V_{\mathit{a}}={\left(C^{\top }C\right)}^{-1}$ and let $\mathit{b}=B\mathit{a}$, where B is $n\times n$ and full rank, so that b can be regarded as an alternative parametrization of the problem. For the same experimental design, the variance matrix associated with the estimates b is given by

$$V_{\mathit{b}}=B{V}_{\mathit{a}}{B}^{\top }.$$

We have

$$|V_{\mathit{b}}|=|B{V}_{\mathit{a}}{B}^{\top }{|=|B|}^2|V_{\mathit{a}}|,$$

so that a strategy that is optimal for the D-measure for parametrization a will also be optimal for parametrization b. For the A-measure,

$$\mathrm{Tr}\left(V_{\mathit{b}}\right)=\mathrm{Tr}\left(B{V}_{\mathit{a}}{B}^{\top }\right)=\mathrm{Tr}\left(B^{\top }B{V}_{\mathit{a}}\right)=\mathrm{Tr}\left(V_{\mathit{a}}{B}^{\top }B\right),$$

and, in general, there is no constant K such that $\mathrm{Tr}\left(B{V}_{\mathit{a}}{B}^{\top }\right)=K\mathrm{Tr}\left(V_{\mathit{a}}\right)$, unless $B^{\top }B$ is a multiple of the identity matrix. In a similar vein, the D-measure is invariant with respect to changes in unit, the A-measure is not.

The eigenvalues ${\lambda }_j$ of $V_{\mathit{a}}$ are related to the size of the confidence ellipses centred about the parameter estimates a. In fact, for a given confidence level, the lengths of the semi-axes are proportional to $\sqrt{{\lambda }_j}$. Hence, the D-measure is related to the square of the volume of the confidence ellipse, while the A-measure is related to the sum of the squares of the semi-axis lengths.

4. Choosing a D-Optimal Subset of n Measurements

Let C be an $m\times n$ observation matrix, $m>n$, representing m measurements that could potentially be made. If we are constrained only to use n measurements, which n should we choose, in order to minimise an aggregate measure of the uncertainties associated with the fitted parameters? It is assumed that there is at least one $n\times n$ submatrix constructed from the rows of C that is full rank, allowing all the parameters a to be estimated. In this section, we consider two algorithms described below in Section 4.1 and Section 4.2, for addressing the issue of subset selection with a view to providing a D-optimal solution. Both algorithms were designed to address other problems but can be repurposed for design of experiment applications.

Let $C_k$ be a full rank $n\times n$ submatrix of C defined by a subset of the rows of C indexed by row indices $I_k=\{i_1,\dots ,i_n\}$, and set $V_k=H_k^{-1}$, where $H=C_k^{\top }{C}_k$. Then,

$$|V_k|=|H_k{|}^{-1}={|C_k|}^{-2},$$

so that choosing $I_k$ to minimise $|V_k|$ amounts to choosing $I_k$ to maximise $|C_k|$. For any square $n\times n$ matrix A with rows ${\mathit{a}}_i$, $|A|$ is the volume of the parallelopiped spanned by the vectors ${\mathit{a}}_i$ in ${\mathbb{R}}^n$. This volume will depend on the Euclidean lengths $\parallel {\mathit{a}}_i{\parallel }_2$ of the vectors and also their mutual independence. For vectors of equal length, $|A|$ is maximised if the vectors ${\mathit{a}}_i$ are orthogonal to each other (and can therefore be rotated or reflected to align with the coordinate axes). In terms of the observation matrix C, the lengths of the row vectors ${\mathit{c}}_i$ reflect the accuracy associated with the measurement $y_i$—more accurate measurements will be accorded a greater weight and therefore larger row vector ${\mathit{c}}_i$ in terms of Euclidean length. The orthogonality of the row vectors relates to how independent the sets of information represented by the measurements are from each other. Consider

$$\left[\begin{array}{c}{y}_1\\ y_2\end{array}\right]=C\left[\begin{array}{c}{a}_1\\ a_2\end{array}\right]+\left[\begin{array}{c}{ϵ}_1\\ {ϵ}_2\end{array}\right],C=C\left(\theta \right)=\left[\begin{array}{cc}1& 0\\ \cos \theta & \sin \theta \end{array}\right]{ϵ}_1,{ϵ}_2\in \mathcal{N}(0,1).$$

Both measurements $y_1$ and $y_2$ are associated with the equal information gain, since $\parallel {\mathit{c}}_1\parallel =\parallel {\mathit{c}}_2\parallel =1$. The data point $y_1$ provides information about $a_1$ alone, while the data point $y_2$ potentially provides information about both $a_1$ and $a_2$. For $\theta \approx 0$, the observation equation involving $y_2$ only serves to provide more information about $a_1$, and $a_2$ remains poorly estimated with $|C(\theta )|=|\sin \theta |\approx 0$. For $\theta \approx \pi /2$, the information provided by $y_2$ provides significant new information about $a_2$ and $|C(\theta )|=|\sin \theta |\approx 1$. In this way, amongst a set of measurements of comparable accuracy, we look for subsets that are mutually maximally independent. Both of the algorithms described below involve trying to select a set of measurements that are sufficiently independent.

4.1. Subset Selection Using the QR Factorisation with Column Pivoting

The QR factorisation of an $m\times n$ matrix $C=QR$ has the following geometrical interpretation. The column vectors ${\mathit{c}}_j$ of C define an n dimensional subspace $\mathcal{C}$ of ${\mathbb{R}}^m$. The orthogonal matrix Q defines an axis system for ${\mathbb{R}}^m$, such that the n columns of $Q_1$ define an axis system for $\mathcal{C}$ and the $m-n$ columns of $Q_2$ define an axis system for the space of vectors ${\mathcal{C}}^{\perp }$ orthogonal to $\mathcal{C}$. The columns for $Q_1$ are constructed so that ${\mathit{q}}_1$ is aligned with ${\mathit{c}}_1$, and so that there is an $r_{11}$ such that ${\mathit{c}}_1=r_{11}{\mathit{q}}_1$. The vector ${\mathit{q}}_2$ is chosen to lie in the plane defined by ${\mathit{c}}_1$ and ${\mathit{c}}_2$, and so there are scalars $r_{12}$ and $r_{22}$ such that ${\mathit{c}}_2=r_{12}{\mathit{q}}_1+r_{22}{\mathit{q}}_2$, etc. This gives the factorisation

$$\left[{\mathit{c}}_1\cdots {\mathit{c}}_n\right]=\left[{\mathit{q}}_1\cdots {\mathit{q}}_n\right]\left[\begin{array}{cccc}{r}_{11}& r_{12}& \cdots & r_{1n}\\ 0& r_{22}& \cdots & r_{2n}\\ 0& 0& \ddots & ⋮\\ 0& 0& \cdots & r_{nn}\end{array}\right],$$

i.e., in matrix notation $C=Q_1{R}_1$.

An equivalent geometrical interpretation is that the orthogonal matrix $Q^{\top }$, a combination of rotations and reflections, transforms C so that its first column is aligned with the first axis in ${\mathbb{R}}^m$, the second column lies in the plane defined by the first and second axes in ${\mathbb{R}}^m$, the third column lies in the three-dimensional space defined by the first three axes in ${\mathbb{R}}^m$, and so on: $Q^{\top }C=R$.

For the standard QR algorithm, the columns of C are processed in the order they appear in the matrix C. For the QR algorithm with column pivoting [23], the order of the columns is changed to improve the numerical properties of the algorithm. The column of largest Euclidean length is aligned with the first axis. Of the remaining columns, the column that is furthest from the space defined by the first axis is then rotated to the plane defined by the first and second axes. At the $k+1$th stage, it is the column furthest from the subspace defined by the first k axes directions that is chosen to be transformed to lie in the space defined by the first $k+1$ axes directions. The output of the algorithm includes the specification of a permutation matrix P that permutes the columns of C, along with an orthogonal Q and upper triangular R such that $CP=QR$.

The first algorithm considered here is based on a subset selection algorithm described in [23] (Section 5.5). The issue addressed by the subset selection algorithm is that of trying to explain a measured response vector y in terms of basis vectors, usually representing covariates x and functions of covariates, stored in an $m\times p$ matrix C. It is assumed that y lies close to a n dimensional subspace defined by a subset of the columns of C. The problem is to estimate n and choose n columns of C that best explain the response, using the QR with column pivoting to select columns that are optimally mutually independent. This information can then be used to predict other, possibly future, responses based on knowledge of the variates represented by the selected columns. The algorithm uses the singular value decomposition (SVD) to estimate n and the QR algorithm with column pivoting to select an appropriate set of n columns.

Our design of experiment problem is different. The observation matrix C of all potential measurements is assumed to be full rank, so estimating the rank of C is not required. More importantly, our selection problem is to select from the rows of C, not the columns. However, the selection of a set of vectors that are maximally mutually independent is directly relevant. Suppose C has QR factorisation $C=Q_1{R}_1$, where $Q_1$ is $m\times n$ and $R_1$ in $n\times n$ and upper-triangular and full rank. The matrix $R_1$ can be used to re-parametrize the problem in terms of $\mathit{b}=R_1^{-1}\mathit{a}$. The observation associated with parameters b is $Q_1$. Since D-optimal designs are independent of parametrization, it is sufficient to work with $Q_1$ as follows. We use column pivoting to find the QR factorisation of $Q_1^{\top }$ so that $Q_1^{\top }P=UT$, where U is an $n\times n$ orthogonal matrix and T is upper-triangular. The first n columns selected by the QR decomposition of $Q_1^{\top }$ specified by P also specify the n rows of C that are a good choice for the D-optimal design. This algorithm is summarised below (Algorithm 1). Its computational complexity is $O\left(m{n}^2\right)$.    

Algorithm 1: SSQR: subset selection algorithm using QR with column pivoting to determine a favourable subset of an observation matrix in terms of the D-measure

Data: An $m\times n$ observation matrix C. Result: Index set $I\subset \{i=1,\dots ,m\}$ where the first n elements specify $C_n$ that represent a good selection of the rows of C in terms of maximising the determinant. Steps: Determine the QR factorization of $C=Q_1{R}_1$, where $Q_1$ is an $m\times n$ orthogonal matrix. Use column pivoting to find the QR factorisation of $Q_1^{\top }$ so that $Q_1^{\top }P=UT$, where U is an $n\times n$ orthogonal matrix and T is upper-triangular. Assign the index I by assigning $I\left(k\right)=i_k$ where $i_k$ specifies the row index of the kth column of P such that $P(i_k,k)=1$.

The SSQR algorithm is not targeted directly at determining a D-optimal design but is an approach that, at least heuristically, has elements that makes it likely to produce a design that is good in terms of the D-measure. The algorithm is an example of a sequential ‘greedy algorithm’, in that at each step, a choice is made that maximises some measure, irrespective of what previous steps have been made, and takes no account of what choices could arise in a future step. Once a row has been selected, there is no opportunity for the row to be deselected at a future iteration. The algorithm is largely independent of the ordering of the rows of the observation matrix C and, other than for reasons of rounding error effects and symmetries in the observation matrix, will arrive at the same solution for different row orderings. As has been mentioned, there is no guarantee that this solution is globally optimal. Example calculations using the algorithm are discussed in Section 5.

4.2. Exchange Approach Based on the Gu and Eisenstat Algorithm

The second algorithm described in this paper is directly targeted at determining a D-optimal design. It is based on the iterative algorithm of Gu and Eisenstat [25] for determining a strongly rank-revealing QR decomposition. Below, we show how it can be adapted to provide a reduction in the D-measure at each iteration through exchanging rows and to stop when no reduction is possible. It starts with a selection $C_0$ of n rows of the $m\times n$ observation matrix C, where $C_0$ is already full rank; e.g., that provided by the SSQR subset selection algorithm of Section 4.1.

We recall that minimising the D-measure is equivalent to maximising the absolute value of the determinant of the upper-triangular factor of the observation matrix. The idea of the algorithm is as follows. Partition $C^{\top }$ as $C^{\top }=\left[\begin{array}{cc}A& B\end{array}\right]$, where A is $n\times n$ and B is $n\times (m-n)$. We assume that $A=C_0^{\top }$ is full rank. The Gu–Eisenstat (GE) algorithm uses the following result:

Proposition 1.

For $n\times n$ full rank matrix A and $n\times p$ matrix B, let $\tilde{A}$ be the matrix constructed by replacing the ith column of A, $1\le i\le n$, by the jth column of B, $1\le j\le p$. Then,

$$|\tilde{A}|=F_{i,j}|A|$$

where $F=A^{-1}B$.

Proof. 

For n-vectors u and v,

$$|A+\mathit{u}{\mathit{v}}^{\top }|=|A|(1+{\mathit{v}}^{\top }{A}^{-1}\mathit{u});$$

(5)

see Appendix A.2. The result follows for $\mathit{v}={\mathit{e}}_i$, the ith column of the identity matrix and $\mathit{u}={\mathit{b}}_j-{\mathit{a}}_i$, noting that $A^{-1}{\mathit{a}}_i={\mathit{e}}_i$.    □

This suggests the following approach, assuming that the first n rows of C are full rank (Algorithm 2).

Algorithm 2: Basic steps to determine a D-optimal subset of an observation matrix

The main computational element of this algorithm is forming $F=A^{-1}B$ at each iteration. Forming F explicitly in this way involves $O\left(n^3\right)$ steps to evaluate $A^{-1}$ and a further $O\left(mn\right)$ steps to evaluate the matrix product. Let Q be orthogonal such $\tilde{A}=Q^{\top }A$ is upper-triangular and set

$$Q^{\top }\left[\begin{array}{cc}A& B\end{array}\right]=\left[\begin{array}{cc}\tilde{A}& \tilde{B}\end{array}\right].$$

(6)

Then, $F=A^{-1}B={\tilde{A}}^{-1}\tilde{B}$, and the latter can be computed efficiently using back-substitution. Thus, step 4 in Algorithm 2 can be replaced by

4.1

Determine orthogonal Q such that $Q^{\top }A$ is upper-triangular.

4.2

Set

$$\left[\begin{array}{cc}A& B\end{array}\right]:=Q^{\top }\left[\begin{array}{cc}A& B\end{array}\right],F=A^{-1}B.$$

Since $|Q^{\top }A|=|A|$, up to sign, the modified algorithm still determines the optimal index set. As described in the Gu and Eisenstat paper, these computations can be made much more efficient by exploiting the fact that at each iteration we are performing a rank-one modification to update an upper triangular matrix in step 2 in Algorithm 2 above. Suppose we have determined to replace ith column of the upper-triangular matrix A by the jth column of B. First, we move the ith column of A to the nth column, as shown schematically in (7).

$$\begin{array}{cccccc}*& *& i& *& *& *\\ & *& i& *& *& *\\ & & i& *& *& *\\ & & & *& *& *\\ & & & & *& *\\ & & & & & *\end{array}⟶\begin{array}{cccccc}*& *& *& *& *& i\\ & *& *& *& *& i\\ & & *& *& *& i\\ & & *& *& *& \\ & & & *& *& \\ & & & & *& \end{array}$$

(7)

This step can be written as

$$\left[\begin{array}{ccc}{A}_1& {\mathit{a}}_i& A_2\end{array}\right]\mapsto \left[\begin{array}{ccc}{A}_1& A_2& {\mathit{a}}_i\end{array}\right]=A{\mathsf{\Pi }}_{i,n}$$

where ${\mathsf{\Pi }}_{i,n}$ is the appropriate permutation matrix. The right-hand matrix is upper-triangular except for the sub-diagonal elements stored in $A_2$. Let Q be the orthogonal matrix where $\tilde{A}=Q^{\top }A{\mathsf{\Pi }}_{i,n}$ is upper-triangular and set $\tilde{B}=Q^{\top }B$. Note that $Q^{\top }$ only acts on rows and columns i to n of $A{\mathsf{\Pi }}_{i,n}$ and rows i to n of B. We have that

$$\tilde{F}={\tilde{A}}^{-1}\tilde{B}={\mathsf{\Pi }}_{i,n}^{\top }F.$$

We can now replace the nth column ${\tilde{\mathit{a}}}_n$ of $\tilde{A}$ by the jth column ${\tilde{\mathit{b}}}_j$ of $\tilde{B}$, forming matrices

$$\overline{A}=\tilde{A}+\mathit{u}{\mathit{e}}_n^{\top },\overline{B}=\tilde{B}-\mathit{u}{\mathit{e}}_j^{\top },\mathit{u}={\tilde{\mathit{b}}}_j-{\tilde{\mathit{a}}}_n.$$

Note that the matrix $\overline{A}$ on the left remains upper-triangular. For the next iteration it remains to calculate $\overline{F}={\overline{A}}^{-1}\overline{B}$. Using the Sherman–Morrison formula [26] (Appendix A.1), we have

$${\overline{A}}^{-1}={\left(\tilde{A}+\mathit{u}{\mathit{e}}_n\right)}^{-1}={\tilde{A}}^{-1}-\frac{\tilde{\mathit{u}}{\tilde{\mathit{v}}}^{\top }}{1+{{\mathit{e}}_n}^{\top }\tilde{\mathit{u}}},\tilde{\mathit{u}}={\tilde{A}}^{-1}\mathit{u},\tilde{\mathit{v}}={\tilde{A}}^{-\top }{\mathit{e}}_n,$$

so that

$$\begin{array}{ccc}\hfill {\overline{A}}^{-1}\overline{B}& =& {\tilde{A}}^{-1}\tilde{B}-{\tilde{A}}^{-1}\mathit{u}{\mathit{e}}_j^{\top }-\frac{\tilde{\mathit{u}}\left({\tilde{\mathit{v}}}^{\top }\tilde{B}\right)}{1+{{\mathit{e}}_n}^{\top }\tilde{\mathit{u}}}+\frac{\tilde{\mathit{u}}\left({\tilde{\mathit{v}}}^{\top }\mathit{u}\right){\mathit{e}}_j^{\top }}{1+{{\mathit{e}}_n}^{\top }\tilde{\mathit{u}}},\hfill \\ & =& \tilde{F}+\tilde{\mathit{u}}{\mathit{w}}^{\top },\hfill \end{array}$$

(8)

where

$$\mathit{w}=\left(\frac{{\tilde{\mathit{v}}}^{\top }\mathit{u}}{1+{\mathit{e}}_n^{\top }\tilde{\mathit{u}}}-1\right){\mathit{e}}_j-\frac{{\tilde{B}}^{\top }\tilde{\mathit{v}}}{1+{{\mathit{e}}_n}^{\top }\tilde{\mathit{u}}}.$$

Thus, $\overline{F}$ is a rank-one update of $\tilde{F}$. The expressions above can be considerably simplified. We note that

$$\tilde{\mathit{u}}=\tilde{F}(:,j)-{\mathit{e}}_n,\frac{{\tilde{\mathit{v}}}^{\top }\mathit{u}}{1+{\mathit{e}}_n^{\top }\tilde{\mathit{u}}}-1=-1/\tilde{F}(n,j),\frac{{\tilde{B}}^{\top }\tilde{\mathit{v}}}{1+{{\mathit{e}}_n}^{\top }\tilde{\mathit{u}}}=-\tilde{F}{(n,:)}^{\top }/\tilde{F}(n,j),$$

so that

$$\mathit{w}=-\frac{1}{{\tilde{F}}_{n,j}}(\tilde{F}{(n,:)}^{\top }+{\mathit{e}}_j).$$

Recalling that $\tilde{F}={\mathsf{\Pi }}_{i,n}^{\top }{A}^{-1}B$, we see that from (8), ${\overline{A}}^{-1}\overline{B}$ can be constructed as a rank-one update of the previously calculated $F=A^{-1}B$. We also note that $\tilde{F}(n,j)=F(i,j)>1$ (since the update is only invoked to increase $|A|$) so that the rank-one update of F is numerically stable. The formation of $\tilde{A}$ and $\overline{A}$ requires $O\left(n^2\right)$ steps and the formation of $\overline{F}$ requires $O\left(mn\right)$ steps. These calculations lead to a computationally efficient and numerically stable version (Algorithm 3) of Algorithm 2:

Algorithm 3: GE: efficient algorithm to determine a D-optimal subset of an observation matrix

The algorithm is a sequential descent algorithm, with each exchange producing a decrease in the determinant of the variance matrix by at least a factor of $1/f^2<1$, ensuring that it stops in a finite number of steps. However, there is no guarantee that the algorithm will determine the globally D-optimal subset, only a locally D-optimal solution. A different ordering of the rows of the observation matrix is quite likely to produce a different solution. Consider the following matrix $C\left(a\right)$:

$$C\left(a\right)=\left[\begin{array}{cccc}\hfill 1.00& \hfill 0.00& \hfill 0.00& \hfill 0.00\\ \hfill 0.00& \hfill 1.00& \hfill 0.00& \hfill 0.00\\ \hfill 0.00& \hfill 0.00& \hfill 1.00& \hfill 0.00\\ \hfill 0.00& \hfill 0.00& \hfill 0.00& \hfill a\\ \hfill 0.50& \hfill 0.50& \hfill 0.50& \hfill 0.50\\ \hfill 0.17& \hfill -0.83& \hfill 0.17& \hfill 0.50\\ \hfill 0.17& \hfill 0.17& \hfill -0.83& \hfill 0.50\\ \hfill -0.83& \hfill 0.17& \hfill 0.17& \hfill 0.50\end{array}\right].$$

(9)

The bottom four rows form an orthogonal matrix Q with $|Q|=1$. However, for $a>0.5$, there is no row interchange that will allow the determinant of the first four rows to be increased, and the GE algorithms stops with the first four rows with determinant a, sub-optimal if $a<1.0$. For $0.5<a<1.0$, the SSQR algorithm, by contrast, selects the bottom four rows, in this case, the globally optimal solution.

4.3. Combined SSQR-GE Algorithm

An effective approach in practice is to combine the SSQR and GE algorithms, using the SSQR algorithm to find an initial selection of columns, and to use the GE interchange algorithm to improve on, if possible, this selection. This combined algorithm has the advantages of being largely independent of the row ordering and will deliver a local minimum based on the SSQR algorithm starting point. Again, there is no guarantee that the solution will be a global optimum, but the examples discussed in Section 5 show that the combined algorithm works well.

4.4. Sequential Approach for Choosing a Set of Additional Observations

The algorithms described in Section 4.1 and Section 4.2 are designed to determine a minimal set of observations that produce D-optimal estimates of the model parameters a. We now consider the case in which estimates of a have already been established. Suppose the (full rank) variance matrix $V=L{L}^{\top }={\left(R^{\top }R\right)}^{-1}$ represents the current knowledge about parameters a and up to m (new) measurements are potentially available, modelled according to

$$y_i\sim \mathcal{N}({\mathit{c}}_i^{\top }\mathit{a},1),i=1,\dots ,m.$$

(As before, we assume that, for each i, $y_i$ and ${\mathit{c}}_i$ have been weighted so that the uncertainty associated with $y_i$ is 1.) Which measurements should we make in order to best improve our knowledge about a? We suppose that resources exist to make p more observations. An exhaustive search of all combinations is only feasible for small p and modest m. The sequential algorithms below add the observations one at a time, with the choice at each stage designed to maximise the information. These are further examples of so-called ‘greedy algorithms’.

Let C be the $m\times n$ observation matrix whose ith row is ${\mathit{c}}_i^{\top }$ and $V_i$ the variance matrix associated with a given that the ith measurement (and only the ith measurement) is made. Then,

$$V_i={(R^{\top }R+{\mathit{c}}_i{\mathit{c}}_i^{\top })}^{-1}=V-\frac{V{\mathit{c}}_i{\left(V{\mathit{c}}_i\right)}^{\top }}{1+{\mathit{c}}_i^{\top }V{\mathit{c}}_i},$$

using the Sherman–Morrison formula [26] (Appendix A.1). Let

$${\mathit{f}}_i=V{\mathit{c}}_i,g_i^2={\mathit{c}}_i^{\top }{\mathit{f}}_i={\mathit{c}}_i^{\top }V{\mathit{c}}_i.$$

(10)

Then,

$$V_i=V-\frac{{\mathit{f}}_i{\mathit{f}}_i^{\top }}{1+g_i^2}=V-{\mathit{u}}_i{\mathit{u}}_i^{\top },{\mathit{u}}_i=\frac{1}{\sqrt{1+g_i^2}}{\mathit{f}}_i.$$

(11)

so that $V_i$ is a symmetric rank-one update of V.

The eigenvalues ${\lambda }_{j,i}$ of $V_i$ can be related to the eigenvalues ${\lambda }_j$ of V through the following result [23] (Section 8.1), [27] (p. 94):

Proposition 2.

Suppose $W=V+\rho \mathit{u}{\mathit{u}}^{\top }$, where V is an $n\times n$ symmetric matrix, u is an n-vector with ${\mathit{u}}^{\top }\mathit{u}=1$ and ρ is a constant. Suppose V and W have eigenvalues

$${\lambda }_1\left(V\right)\ge {\lambda }_1\left(V\right)\ge {\lambda }_2\left(V\right)\ge \cdots \ge {\lambda }_n\left(V\right),$$

and

$${\lambda }_1\left(W\right)\ge {\lambda }_1\left(W\right)\ge {\lambda }_2\left(W\right)\ge \cdots \ge {\lambda }_n\left(W\right).$$

If $\rho \ge 0$, then

$${\lambda }_i\left(W\right)\in [{\lambda }_i\left(V\right),{\lambda }_{i-1}\left(V\right)],i=2,\dots ,n,$$

while if $\rho \le 0$, then

$${\lambda }_i\left(W\right)\in [{\lambda }_{i+1}\left(V\right),{\lambda }_i\left(V\right)],i=1,\dots ,n-1.$$

In either case, there exist n constants ${\gamma }_i$, ${\gamma }_i\ge 0$, with ${\gamma }_1+\cdots {\gamma }_n=1$, such that

$${\lambda }_i\left(W\right)={\lambda }_i\left(V\right)+{\gamma }_i\rho .$$

4.4.1. Sequential D-Optimality Algorithm

We use the relationship (5) to write

$$\begin{array}{ccc}\hfill |V_i|& =& |(R^{\top }R+{\mathit{c}}_i{\mathit{c}}_i^{\top }{|}^{-1},\hfill \\ & =& |R^{\top }R{|}^{-1}|{\left(1+{\mathit{c}}_i^{\top }{\left(R^{\top }R\right)}^{-1}{\mathit{c}}_i\right)}^{-1},\hfill \\ & =& |V|{\left(1+{{\mathit{c}}_i}^{\top }V{\mathit{c}}_i\right)}^{-1},\hfill \\ & =& \frac{1}{1+g_i^2}|V|,g_i^2={\mathit{c}}_i^{\top }V{\mathit{c}}_i.\hfill \end{array}$$

A sequential algorithm can therefore be designed which at each stage chooses the i that maximises $g_i^2$. The algorithm can be run in two modes: the first in which the algorithm can choose any row at most one time, the second in which there is no limit in the number of times any row can be selected. Having a row repeated k times is equivalent to having a single more accurate observation with its uncertainty reduced by a factor of $1/\surd k$. In the second mode, it is often the case that additional measurements correspond to repeat measurements chosen from the D-optimal set of n points. In the algorithm below (Algorithm 4), the mode is controlled by the flag $i_R$:

Algorithm 4: Sequential approach to optimising the D-measure in choosing p observations from a possible $m>p$ observations

For each q, the only operations required are matrix–vector multiplications and similar, requiring at most $O\left(mn\right)$ steps. This compares with $O\left(m{n}^3\right)$ steps for an algorithm not exploiting rank-one updates.

4.4.2. Expected Information Gain Calculations

The quantities $t_q$ calculated at each stage give the reduction in determinant of the variance matrix at the qth step. The following argument can be used to provide an expected value for $t_q$, under the assumption that the rows of C each provide approximately the same amount of information. Suppose C consists of rows with a 1 in one column and zeros elsewhere. A D-optimal design would select rows so that the ones are distributed equally amongst the columns. If a total of q rows are selected, the variance matrix $V_q$ is such that

$$V_q\approx \frac{n}{q}I,|V_q|\approx {\left(\frac{n}{q}\right)}^n.$$

From this argument, the expected value of $t_q$, $q>n$, is such that

$$t_q\approx \tilde{t}\left(q\right)={\left(\frac{q-1}{q}\right)}^n.$$

(12)

In terms of the geometric mean $\overline{d}$ of the eigenvalues of the variance matrices, we have

$$\overline{d}\left(V_q\right)=\frac{q-1}{q}\overline{d}\left(V_{q-1}\right).$$

For $q=n+1$, we have

$$\frac{1}{2}\ge {\left(\frac{n}{n+1}\right)}^n\ge \frac{1}{e}\approx 0.37,$$

where the left hand value corresponds to the case $n=1$, and $1/e$ is the limiting value as $n\to \infty$.

4.4.3. Sequential A-Optimality Algorithm

We can also construct a counterpart of Algorithm 4 for A-optimality. For any v, $\mathrm{Tr}\left(\mathit{v}{\mathit{v}}^{\top }\right)={\mathit{v}}^{\top }\mathit{v}={\parallel \mathit{v}\parallel }^2$. With ${\mathit{f}}_i$ and $g_i$ as in (10), if

$${\tau }_i^2=f_i^2/(1+g_i^2),f_i^2={\mathit{f}}_i^{\top }{\mathit{f}}_i,$$

then

$$\mathrm{Tr}\left(V_i\right)=\mathrm{Tr}\left(V\right)-{\tau }_i^2.$$

This last expression calculates the change in the A-measure for a given choice of one additional measurement. A sequential algorithm can therefore be designed which at each stage chooses the i that maximises ${\tau }_i^2$. Each stage requires the calculation of the quantities $f_i^2$ and $g_i^2$, $i=1,\dots ,m$. If V is the current variance matrix, then the $f_i^2$ are given by the sum of the squares of row elements of $F=V{C}^{\top }$. Suppose measurement k maximises ${\tau }_i^2$. Then, the update of F is given by

$$F=(V-{\mathit{u}}_k{\mathit{u}}_k^{\top })C^{\top }=F-{\mathit{u}}_k{\mathit{w}}^{\top },\mathit{w}=C{\mathit{u}}_k,$$

and the updated of $g_i^2$ is given by $g_i^2-w_i^2$, with the latter evaluated as $(g_i-w_i)(g_i+w_i)$ for better numerical accuracy. These steps are summarised in (Algorithm 5):

Algorithm 5: Sequential approach to optimising the A-measure in choosing p observations from a possible $m>p$ observations

As for Algorithm 4, the computational requirement for each step is $O\left(mn\right)$.

5. Numerical Examples

5.1. Polynomial Calibration Curves

Many instrument response functions are modelled as polynomial curves, most commonly, of order 2 (degree 1), i.e., modelling a linear response, but higher orders arise, e.g., in the calibration of platinum resistance thermometers [28], or in the calibration of a stage motion [29] where polynomials of order 10 are involved. (Design problems for problems for polynomial models augmented by Gaussian process models are considered in [4].) Suppose a polynomial response function of order n is defined on the interval [−1,1], and that it is desired to calibrate the response function using n calibration points $x_i\in [-1,1]$. Which choice of n points is D-optimal? For a linear (quadratic) response $n=2$ ($n=3$), it is intuitively clear that the end points $x=\pm 1$ (and $x=0$) are optimal. For higher orders, the optimal choice is not obvious. In fact, the D-optimal set [30] for order n is given by solutions of

$$(1-x^2){\dot{L}}_n\left(x\right)=0,$$

(13)

where ${\dot{L}}_n\left(x\right)$ is the derivative of Legendre polynomial of order n (degree $n-1$). These solutions can be approximated by the so-called ‘arcsine’ points

$$x_i=\cos \left(\pi \frac{n-1-i}{n-1}\right),i=0,1,\dots ,n-1.$$

(14)

Table 1. D-optimal calibration points for polynomial calibration curves of order $n=4,\dots ,11$ defined in the interval $[-1,1]$. For each set of four columns, the first column, labelled $x_0$, gives the arcsine estimates of the calibration points in (14); the second column, labelled $x_{\mathrm{QR}}$, gives the optimal the calibration points determined by the SSQR algorithm; the third column, labelled $x_{\mathrm{GE}}$, gives the estimates determined by the combined SSQR-GE algorithm; while the fourth column, labelled $x^{*}$, gives the optimal the calibration points derived from (13).

${\mathit{x}}_{\mathbf{0}}$	${\mathit{x}}_{\mathrm{QR}}$	${\mathit{x}}_{\mathrm{GE}}$	${\mathit{x}}^{*}$	${\mathit{x}}_{\mathbf{0}}$	${\mathit{x}}_{\mathrm{QR}}$	${\mathit{x}}_{\mathrm{GE}}$	${\mathit{x}}^{*}$
$n=4$				$n=5$
−1.000	−1.000	−1.000	−1.000	−1.000	−1.000	−1.000	−1.000
−0.500	−0.488	−0.447	−0.447	−0.707	−0.669	−0.655	−0.655
0.500	0.437	0.447	0.447	0.000	0.006	0.000	0.000
1.000	1.000	1.000	1.000	0.707	0.686	0.655	0.655
				1.000	1.000	1.000	0.000
$n=6$				$n=7$
−1.000	−1.000	−1.000	−1.000	−1.000	−1.000	−1.000	−1.000
−0.809	−0.786	−0.765	−0.765	−0.866	−0.845	−0.830	−0.830
−0.309	−0.286	−0.285	−0.285	−0.500	−0.484	−0.469	−0.469
0.309	0.308	0.285	0.285	0.000	0.002	0.000	0.000
0.809	0.779	0.765	0.765	0.500	0.493	0.469	0.469
1.000	1.000	1.000	1.000	0.866	0.841	0.830	0.830
				1.000	1.000	1.000	0.000
$n=8$				$n=9$
−1.000	−1.000	−1.000	−1.000	−1.000	−1.000	−1.000	−1.000
−0.901	−0.880	−0.872	−0.872	−0.924	−0.908	−0.900	−0.900
−0.623	−0.613	−0.592	−0.592	−0.707	−0.692	−0.677	−0.677
−0.223	−0.211	−0.209	−0.209	−0.383	−0.383	−0.363	−0.363
0.223	0.225	0.210	0.209	0.000	−0.002	0.000	0.000
0.623	0.608	0.592	0.592	0.383	0.376	0.363	0.363
0.901	0.882	0.872	0.872	0.707	0.695	0.677	0.677
1.000	1.000	1.000	1.000	0.924	0.906	0.900	0.900
				1.000	1.000	1.000	0.000
$n=10$				$n=11$
−1.000	−1.000	−1.000	−1.000	−1.000	−1.000	−1.000	−1.000
−0.940	−0.925	−0.920	−0.920	−0.951	−0.938	−0.934	−0.934
−0.766	−0.753	−0.739	−0.739	−0.809	−0.796	−0.784	−0.784
−0.500	−0.493	−0.478	−0.478	−0.588	−0.580	−0.565	−0.565
−0.174	−0.177	−0.165	−0.165	−0.309	−0.311	−0.296	−0.296
0.174	0.168	0.165	0.165	0.000	−0.001	0.000	0.000
0.500	0.497	0.478	0.478	0.309	0.307	0.296	0.296
0.766	0.751	0.739	0.739	0.588	0.582	0.566	0.565
0.940	0.925	0.920	0.920	0.809	0.795	0.785	0.784
1.000	1.000	1.000	1.000	0.951	0.939	0.934	0.934
				1.000	1.000	1.000	0.000

gives the results of applying the SSQR algorithm and the combined SSQR-GE algorithm to estimate the D-optimal calibration points for the cases $n=4,5,\dots ,11$ choosing from $m=2001$ possible points equally spaced in the interval $[-1,1]$ in steps of 0.001. For these calculations, we have used Chebyshev polynomial basis functions [31] for numerical stability (but the optimal choice of points is independent of the choice of basis functions.) The arcsine estimates are also given in the table, along with the D-optimal points derived from (13); see Appendix B. The results show that, in terms of closeness to the D-optimal points, the SSQR algorithm improves on the arcsine estimates, while the GE algorithm improves the SSQR solution and is accurate to three decimals in almost all cases (there is one exception for the case $n=11$).

Table 2. For $n=4,5,\dots ,11$, rows two to five give the computed geometric mean measure $\overline{d}\left(V_{\mathit{a}}\right)$ given in (4) for evenly spaced calibration points, the arcsine solutions, the SSQR solutions, and the SSQR-GE solutions given in Table 1. The final row gives the number of GE exchanges undertaken to improve the SSQR solutions.

n	4	5	6	7	8	9	10	11
${\overline{d}}_0$	0.4871	0.4152	0.3748	0.3511	0.3379	0.3316	0.3304	0.3332
${\overline{d}}_{\mathrm{AS}}$	0.4714	0.3789	0.3175	0.2734	0.2403	0.2143	0.1935	0.1763
${\overline{d}}_{\mathrm{QR}}$	0.4682	0.3746	0.3130	0.2691	0.2362	0.2107	0.1901	0.1733
${\overline{d}}_{\mathrm{GE}}$	0.4673	0.3735	0.3119	0.2682	0.2354	0.2099	0.1894	0.1726
$n_{\mathrm{GE}}$	3	6	7	9	12	17	19	21

gives the computed geometric mean measure $\overline{d}\left(V_{\mathit{a}}\right)$ given in (4) for evenly spaced calibration points, the arcsine solutions, the SSQR solutions, and the SSQR-GE solutions given in Table 1. Also shown is the number of GE exchanges undertaken to improve the SSQR solutions. The table reflects the improvements provided by the SSQR and SSQR-GE algorithms. The number of GE exchanges is small compared to $m=2001$. The evenly spaced calibration points are markedly worse than the other solutions for larger n.

Figure 1 graphs the values of $t_q$, $q=1,\dots ,100$, and its analytical approximation $\tilde{t}\left(q\right)$ given by (12), determined by the sequential D-optimal Algorithm 4, without repeat measurements, for the cases $n=4$, upper curve, and $n=11$, lower curve, starting with the SSQR-GE solution. The values of $t_q$ are batched with a significant jump after every n additional measurement. Without the repeat measurement, the algorithm tends to select available points closest to the D-optimal n points. If the algorithm is run with repeat measurements allowed, the algorithm selects the same set of D-optimal n points again and again. This behaviour is reflected in the graphs of $t_q$ for this case; see Figure 2.

5.2. Tensor Product of Polynomials

Our second example arose from a problem relating to the estimation of heat distribution on an plate [32]. The application involved modelling the heat distribution as a Gaussian process model [33] involving tensor product polynomials, along with spatially correlated effects. A sub-problem is related to choosing a set of calibration points from a $10\times 14$ grid of possible sensor locations in order to determine the tensor product mean function. The tensor product model is of the form

$$f(x,y,\mathit{a})=\sum _{i=1}^{n_x}\sum _{j=1}^{n_y}{a}_{ij}{f}_i\left(x\right)g_j\left(y\right),$$

where $f_i\left(x\right)$ and $g_j\left(y\right)$ are Chebyshev polynomial basis functions. Figure 3 shows the solution calibration points for $n_x=n_y=5$ constructed from polynomial basis functions up to degree 4. Two solutions are shown: the first in which the calibration points are selected from a $131\times 91$ finer grid of 11,921 points, the second in which the calibration points are selected from a $14\times 10$ grid of 140 points. The solutions found for the fine grid coincide with a grid $5\times 5$ grid of points $(x_k^{*},y_l^{*})$, where the $x_k^{*}$, $k=1,\dots ,5$ are the D-optimal solutions derived from the solution of (13) on the interval [0, 20] and the $y_l^{*}$ are the D-optimal solutions on the interval [0, 10]. For this arrangement, the observation matrix C associated with the polynomial fit is the tensor product $C=C_x\otimes C_y$. For square, full rank matrices A and B of size $n_A\times n_A$ and $n_B\times n_B$

$$|A\otimes B|={|A|}^{n_A}{|B|}^{n_B}.$$

Hence, it is consistent that the D-optimal points for the tensor product are constructed from the D-optimal points for the x- and y-basis functions. The solution points for the coarse grid are seen to approximate the solution for the fine grid, given the constraints imposed by the coarse grid. The coarse grid solution points do not represent a regular grid in this case. For the case of the fine grid, the GE algorithm took 16 iterations to reach a local optimum, starting with the SSQR solution. For the coarse grid, the GE algorithm confirmed that the SSQR solution is locally optimal.

Figure 3. Tensor product polynomial calibration points found by the SSQR-GE algorithm. The points marked with a ‘+’ are those selected from a fine grid, those marked ‘o’ from a coarse grid.

Figure 3. Tensor product polynomial calibration points found by the SSQR-GE algorithm. The points marked with a ‘+’ are those selected from a fine grid, those marked ‘o’ from a coarse grid.

Figure 4 and Figure 5 plot the values of $t_q$, $q=1,\dots ,200$, and its analytical approximation $\tilde{t}\left(q\right)$ given by (12), determined by the sequential D-optimal Algorithm 4, for the fine and coarse grids, respectively, with and without repeat measurements. For a coarse grid and no repeat measurements, the additional calibration points that are permitted consistently provide less information compared to the expected information gain given by $\tilde{t}$ in (12). For a coarse grid of just 140 points, later additional measurements have to be selected from locations that are less informative in terms of improving the D-measure. With repeat measurements on the coarse grid, the actual information gain is in line with the expected information gain.

5.3. Coordinate Metrology

In coordinate metrology, the size and shape of a geometric surface $\mathcal{S}$ is determined from the measurement of the coordinates $x_i$ of points lying on the surface. If the geometric surface $\mathcal{S}=\mathcal{S}\left(\mathit{a}\right)$ is parametrized by $\mathit{a}={(a_1,\dots ,a_n)}^{\top }$ and $d(\mathit{x},\mathit{a})$ is a measure of the distance from a point x to $\mathcal{S}\left(\mathit{a}\right)$, then for measurement points ${\mathit{x}}_{1:m}=\{{\mathit{x}}_i,i=1,\dots ,m\}$, $m\ge n$, estimates of a are determined by solving the least squares optimisation problem

$$\underset{\mathit{a}}{min}=\sum _{i=1}^m{d}^2({\mathit{x}}_i,\mathit{a}).$$

If the variance matrix associated with the measured points coordinates is ${\sigma }^{2}I$, then the variance matrix $V_{\mathit{a}}$ associated with the fitted parameters is approximated by

$$V_{\mathit{a}}={\sigma }^2{\left(J^{\top }J\right)}^{-1},J_{ij}=\frac{\partial d}{\partial a_j}({\mathit{x}}_i,\mathit{a}),$$

involving the Jacobian matrix J. Thus, it is possible to determine optimal measurement strategies, i.e., where to measure on the surface $\mathcal{S}$, in terms of the determinant of the Jacobian matrix J. Often, the case of choosing a minimal set of points with $m=n$ is of interest. For simple geometries, an optimal set of points is straightforward to arrive at, e.g., three points equally spaced around a circle, but for more complex geometries, the choice of an optimal set may not be obvious, as considered below.

Figure 6 shows a reference artefact designed and calibrated by the Czech Metrology Institute (CMI) [34,35]. The form of the geometric surface is a hyperbolic paraboloid given by the equation

$$z=24+(x/8-1)(y/8-1).$$

(15)

A paraboloid with its axis approximately parallel to the z-axis is parametrized in terms of two rotation angles $\alpha$ and $\beta$ and six further parameters b such that

$$\widehat{z}=b_1{\widehat{x}}^2+b_2{\widehat{y}}^2+b_3\widehat{x}\widehat{y}+b_4\widehat{x}+b_5\widehat{y}+b_6,\widehat{\mathbf{x}}=R(\alpha ,\beta )\mathbf{x},$$

where $R(\alpha ,\beta )=R_x\left(\alpha \right)R_y\left(\beta \right)$ is the rotation matrix constructed from the plane rotation matrices

$$R_x\left(\alpha \right)=\left[\begin{array}{ccc}\hfill 1& \hfill 0& \hfill 0\\ \hfill 0& \hfill \cos \alpha & \hfill -\sin \alpha \\ \hfill 0& \hfill \sin \alpha & \hfill \cos \alpha \end{array}\right],R_y\left(\beta \right)=\left[\begin{array}{ccc}\hfill \cos \beta & \hfill 0& \hfill \sin \beta \\ \hfill 0& \hfill 1& \hfill 0\\ \hfill -\sin \beta & \hfill 0& \hfill \cos \beta \end{array}\right].$$

For the surface (15), nominally ${\mathit{b}}^{\top }=(0,0,1/64,-1/8,-1/8,25)$. Figure 7, Figure 8 and Figure 9 show the distribution of the x- and y-coordinates for three sets $X_k$, $k=1,2,3$, of eight calibration points. The first two are given by ‘expert judgement’; that is, possible choices for a measurement strategy based on engineering practice. The first expert choice, $X_1$, Figure 7, represents an approximately uniform distribution, a first guess at a good distribution. This distribution of points is sufficient to determine all the parameters, and the associated $8\times 8$ Jacobian matrix is full rank. (If the x- and y- coordinates of the data points are rotated through 45 degrees about the z-axis, the associated Jacobian matrix is rank deficient by three.) The second set $X_2$ came about through experimenting with moving the inner four points in dataset $X_1$ closer or further from the origin of the $xy$-plane. It turned out that moving the four points further from the origin was better, resulting in dataset $X_2$, Figure 8. At first glance, the distribution of points in $X_2$ does not look a good choice. The third dataset is that obtained using the SSQR-GE algorithm. The GE algorithm performed 11 interchanges starting from the SSQR solution. The optimal choice has seven points along the boundary in the $xy$-plane and one point in the interior. Any interior point not too far from the origin would also result in a good selection. In terms of the aggregate measure ${\overline{d}}_k$ defined in (4) associated with the three datasets, we have ${\overline{d}}_1=6.0{\overline{d}}_3$ and ${\overline{d}}_2=2.2{\overline{d}}_3$.

Table 3. Standard uncertainty factors $u\left(a_j\right)/\sigma$ associated with the hyperbolic paraboloid parameters $\mathit{a}={(\alpha ,\beta ,{\mathit{b}}^{\top })}^{\top }$ for three eight-point datasets $X_k$, $k=1,2,3$.

$\mathit{u}\left({\mathit{a}}_{\mathit{j}}\right)/\mathit{\sigma }$	$\mathit{\alpha }$	$\mathit{\beta }$	${\mathit{b}}_1$	${\mathit{b}}_2$	${\mathit{b}}_3$	${\mathit{b}}_4$	${\mathit{b}}_5$	${\mathit{b}}_6$
$X_1$	2.3867	2.3867	0.0049	0.0049	0.0139	2.6219	2.6219	10.5328
$X_2$	0.7259	0.7259	0.0042	0.0042	0.0035	0.4999	0.4999	2.8188
$X_3$	0.2049	0.2413	0.0012	0.0014	0.0014	0.2883	0.2520	1.4690

gives the square roots ${\mathit{v}}_i$ of the diagonal elements of ${\left(J_k^{\top }{J}_k\right)}^{-1}$, where $J_k$ is the Jacobian matrix associated with the kth dataset. If the variance matrix associated with the measurements is ${\sigma }^{2}I$, then $v_i=u\left(a_j\right)/\sigma$, where $u\left(a_j\right)$ is the estimate of the standard uncertainty associated with the jth parameter [36,37]. In terms of standard uncertainties, the D-optimal dataset $X_3$ is far better for estimating the surface parameters than the other two datasets.

Figure 10 plots the values of $t_q$, $q=1,\dots ,100$, and its analytical approximation $\tilde{t}\left(q\right)$ given by (12), determined by the sequential D-optimal Algorithm 4, with and without repeat measurements. For higher values of q, there is good agreement between $t_q$ and ${\tilde{t}}_q$.

5.4. Calibration of a Network of Standards Using a Comparator

In this section, we consider the calibration of a number of standards, starting with one calibrated standard and using a comparator to calibrate the other standards relative to the calibrated standard. An example application is in mass calibration using a mass balance [38,39,40]. Given a calibrated standard $A_1$, nominally of 1 kg, and two uncalibrated masses $A_2$ and $A_3$ each of nominal mass 0.5 kg, $A_2$ and $A_3$ can be calibrated using two measurements, the first comparing $A_1$ with the combined mass of $A_2$ and $A_3$, and the second comparing $A_2$ and $A_3$. A second example from length metrology relates to the calibration of gauge blocks [41], where it is possible to align two gauge blocks end to end, a process called ‘wringing’, to define a combined length. This discussion is relevant to the calibration of any set of artefacts defining an extensive quantity where the artefacts can be grouped to define a quantity that is (nominally) the sum of the individual quantities. We use the term ‘network’ to reflect the fact that the estimates of all the quantities are statistically correlated and an experiment that updates the information about one artefact will also update information about the other artefacts.

Let $\mathit{a}={(a_1,\dots ,a_n)}^{\top }$ represent the values of the quantities associated with the artefacts where the nominal values ${\mathit{a}}^{♯}={(a_1^{♯},\dots ,a_n^{♯})}^{\top }$ of the quantities are such that the artefacts can be grouped into sets of two groups, with each group being associated with the same nominal value. Thus, there are disjoint subsets $L_i$, $R_i\subset I=\{1,2,\dots ,n\}$, $L_i\cap R_i=\varnothing$, $i=1,\dots ,m$, such that

$$\sum _{q\in L_i}{a}_q\approx \sum _{q\in R_i}{a}_q.$$

An example considered further below involves a set of artefacts having nominal values 1.00, 0.50, 0.50, 0.20, 0.20, 0.10, 0.10, 0.05, and 0.05. The fact that the two subsets are associated with the same nominal value allows a comparator measurement to be made, with the associated observation equation

$$0\approx y_i=\sum _{q\in L_i}{a}_q-\sum _{q\in R_i}{a}_q+{ϵ}_i,{ϵ}_i\in \mathcal{N}(0,{\sigma }_i^2),$$

where ${ϵ}_i$ represents a random effect associated with the comparator measurement. The design problem is to define measurement strategies encoded by $L_i$ and $R_i$, $i=1,\dots ,n$, that provide the most accurate calibration of the standards. For the example considered here, there are almost 400 viable experiments.

Table 4. Design of a calibration experiment given by ‘expert judgement’.

1.00	0.50	0.50	0.20	0.20	1.0	1.0	0.05	0.05
1	0	0	0	0	0	0	0	0
1	−1	−1	0	0	0	0	0	0
0	1	−1	0	0	0	0	0	0
0	1	0	−1	−1	−1	0	0	0
0	0	1	−1	−1	0	−1	0	0
0	0	0	1	−1	0	0	0	0
0	0	0	1	0	0	0	−1	−1
0	0	0	0	0	1	0	−1	−1
0	0	0	0	0	0	0	1	−1

and

Table 5. Designs of four calibration experiments calculated using the SSQR-GE algorithm with uncertainties ${\sigma }_i$ calculated according to (16) for the different values of ${\sigma }_R$, ${\sigma }_N$, and ${\sigma }_V$ given in Table 6.

1.00	0.50	0.50	0.20	0.20	1.0	1.0	0.05	0.05
1	0	0	0	0	0	0	0	0
1	0	−1	−1	−1	−1	−1	1	1
1	−1	−1	1	−1	1	−1	1	−1
1	−1	−1	1	−1	−1	1	−1	1
1	−1	−1	−1	1	1	−1	−1	1
1	−1	−1	−1	1	−1	1	1	−1
0	1	0	−1	−1	0	0	−1	−1
0	1	−1	0	1	−1	0	−1	−1
0	0	1	−1	0	−1	−1	−1	−1
1	0	0	0	0	0	0	0	0
1	−1	−1	0	0	0	0	0	0
0	1	0	−1	−1	−1	−1	1	1
0	1	−1	1	−1	1	−1	−1	1
0	1	−1	0	1	−1	0	−1	−1
0	1	−1	0	−1	1	1	1	−1
0	0	1	−1	−1	0	0	−1	−1
0	0	0	1	0	−1	−1	1	−1
0	0	0	1	−1	−1	1	−1	1
1	0	0	0	0	0	0	0	0
1	−1	−1	0	0	0	0	0	0
0	1	0	−1	−1	0	−1	0	0
0	1	−1	0	0	0	0	0	0
0	0	0	1	0	−1	−1	0	0
0	0	0	1	−1	0	0	0	0
0	0	0	0	0	1	0	−1	−1
0	0	0	0	0	1	−1	0	0
0	0	0	0	0	0	0	1	−1
1	0	0	0	0	0	0	0	0
1	−1	−1	0	0	0	0	0	0
0	1	0	−1	−1	0	−1	0	0
0	1	−1	0	0	0	0	0	0
0	0	0	1	0	−1	−1	0	0
0	0	0	1	−1	0	0	0	0
0	0	0	0	0	1	0	−1	−1
0	0	0	0	0	1	−1	0	0
0	0	0	0	0	0	0	1	−1

show five designs for the calibration of the $n=9$ standards. A ‘1’ in column j indicates that $j\in L_i$, while a ‘$-1$’ indicates that $j\in R_i$. Each array of $9\times 9$ elements is in fact the unweighted observation matrix C associated with the particular design. The first row in each design represents the observation associated with the calibration of the first standard using an absolute measurement system. This first observation is the only source of information about the first standard and must be included in the optimal solution. The uncertainties associated with the remaining standards scale with the uncertainty associated with the first observation.

The first design, given in Table 4, is one defined by ‘expert judgement’, and represents a typical set of experiments that would be performed in practice. (In fact, for this set of standards, it is not entirely obvious how to choose a design that will calibrate all of the standards). The next design, indicated by the first nine rows in Table 5, is that determined by the SSQR-GE algorithm with ${\sigma }_1={\sigma }_C=1.0$ and ${\sigma }_i={\sigma }_R=0.5$, $i=2,\dots ,n$. The most striking and, at first sight, surprising feature of this second design is that almost all the artefacts are involved in experiments 2 to n. The reason for this is that, in some sense, the comparator uncertainties represented by ${\sigma }_i={\sigma }_R$ are partitioned amongst the artefacts involved, and more artefacts being included in an experiment means that each artefact is assigned a smaller uncertainty.

The two designs considered so far are based on the assumption that all the comparator measurements are associated with the same uncertainty ${\sigma }_R$, irrespective of the number of artefacts and the nominal values of the artefacts. A more plausible assignment of uncertainties would take into account both the number and nominal values. For a mass balance, it would be usual for the associated uncertainty to have a component that varied in proportion to the mass. Similarly, a length measuring device will usually have one or more influencing factors, e.g., refractive index effects for the case of a laser interferometric comparator [41], that vary in proportion to length. The use of multiple artefacts in each experiment will also likely introduce effects that will add to the uncertainties.

For each experiment, let $n_i$ be the total number of artefacts involved and let

$$v_i=\sum _{q\in L_i\cup R_i}{a}_q^{♯},$$

be a measure of the total value of the quantities involved in the ith experiment. Given ${\sigma }_R$, ${\sigma }_V$ and ${\sigma }_N$, we can assign ${\sigma }_i$ according to

$${\sigma }_i^2={\sigma }_R^2+n_{i,2}{\sigma }_N^2+v_i^2{\sigma }_V^2,n_{i,2}=max(n_i-2,0),i=2,\dots ,n.$$

(16)

Table 5 presents the optimal designs calculated using the SSQR-GE algorithm with ${\sigma }_i$ calculated as in (16) for the four different values of ${\sigma }_R$, ${\sigma }_N$ and ${\sigma }_V$ shown in

Table 6. Uncertainties associated with $a_i$ for four different assignments of ${\sigma }_i$, $i=2,\dots ,n=9$. The columns labelled E give the uncertainties associated with ‘expert judgement’ design given in Table 4, while the columns labelled ‘O’ give those associated with the optimal designs, given, in Table 5, calculated by the SSQR-GE algorithm. Also shown in the table are the aggregate measure $\overline{d}={|V_{\mathit{a}}|}^{1/n}$ of uncertainty, the total number ${\sum }_i{n}_i$ of artefacts, and a measure ${\sum }_i{v}_i$ of the total nominal values involved in each set of experiments. The last three rows give the values of ${\sigma }_R$, ${\sigma }_N$, and ${\sigma }_V$ used to calculate ${\sigma }_i$, $i=2,\dots ,n$.

${\mathit{a}}_{\mathit{i}}^{♯}$	E	O1	E	02	E	O3	E	04
1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
0.50	0.61	0.56	0.66	0.64	0.69	0.69	1.04	1.04
0.50	0.61	0.55	0.66	0.64	0.69	0.69	1.04	1.04
0.20	0.39	0.31	0.43	0.40	0.60	0.58	0.50	0.57
0.20	0.49	0.29	0.52	0.39	0.61	0.58	0.54	0.60
0.10	0.57	0.26	0.61	0.32	0.90	0.44	0.57	0.36
0.10	0.91	0.27	1.03	0.36	1.64	0.45	1.34	0.34
0.05	0.35	0.20	0.36	0.28	0.40	0.48	0.29	0.27
0.05	0.35	0.20	0.36	0.28	0.40	0.48	0.29	0.27
$\overline{d}$	0.17	0.06	0.21	0.12	0.21	0.13	0.21	0.15
$\sum v_i$	7.0	17.4	7.0	11.0	7.0	6.3	7.0	6.3
$\sum n_i$	24	62	24	48	24	22	24	22
${\sigma }_R$	0.50	0.50	0.50	0.50	0.20	0.20	0.20	0.20
${\sigma }_N$	0.00	0.00	0.20	0.20	0.80	0.80	0.20	0.20
${\sigma }_V$	0.00	0.00	0.20	0.20	0.20	0.20	0.80	0.80

. The third design in Table 5 is the result penalising experiments with a large number of artefacts, while the fourth design in the table is the result of penalising experiments with a larger value of $v_i$.

Table 6 gives the uncertainties $u\left(a_i\right)$ associated with $a_i$ for four different assignments of ${\sigma }_i$, $i=2,\dots ,n=9$ for the designs in Table 4 and Table 5. Also shown in the table is the aggregate measure $\overline{d}={|V_{\mathit{a}}|}^{1/n}$ of uncertainty, the total number ${\sum }_i{n}_i$ of artefacts and a measure ${\sum }_i{v}_i$ of the total nominal values involved in each set of experiments. The last three rows give the values of ${\sigma }_R$, ${\sigma }_N$ and ${\sigma }_V$ used to calculate ${\sigma }_i$, $i=2,\dots ,n$. In all cases, the SSQR-GE algorithm leads to significant improvements over the expert judgement design.

6. Concluding Remarks

In this paper, we have looked at the experimental design problem: given m potential observations to determine n parameters, $m>n$, what is the best choice of n observations. In the context of least squares estimation, the problem was formulated as finding the $n\times n$ submatrix of the complete $m\times n$ observation matrix that has maximum determinant, corresponding to the D-optimality criterion. We described two algorithms, the SSQR and GE algorithms, to address this problem. Both were adapted from numerical linear algebra algorithms associated with the QR factorisation of a matrix. We also described algorithms for updating estimates of the problem parameters that are sequentially optimal with the respect to the D-optimality and A-optimality criteria. We illustrated the behaviour of the algorithms on a number of applications drawn from the field of metrology. These algorithms are straightforward to implement. The SSQR algorithm, in particular, is a minor adaption of the standard QR factorisation algorithm with pivoting. The algorithms enable computationally efficient implementations, exploiting rank-one matrix updating techniques. The algorithms are descent-type algorithms that will converge to a local minimum that is not necessarily a global minimum. In the examples considered, the algorithms determined effective experimental designs, often with much better performance than ‘hand-crafted’ designs based on expert judgement.