The Polytope of Optimal Approximate Designs: Extending the Selection of Informative Experiments

Let $\mathcal{X}\subset\mathbb{R}^{k}$ be a compact design space, that is, the set of all design points available for the trials. Note that such a definition of $\mathcal{X}$ includes both continuous and discrete spaces.

Consider a continuous function $\mathbf{f}:\mathcal{X}\to\mathbb{R}^{m}$, $m\geq 1$, that defines a linear regression model on $\mathcal{X}$ in the following sense: provided that the experimental conditions of a trial are determined by the design point $\mathbf{x}\in\mathcal{X}$, the observation is $Y_{\mathbf{x}}=\mathbf{f}^{\prime}(\mathbf{x})\beta+\epsilon$, where $\beta\in\mathbb{R}^{m}$ is a vector of unknown parameters, and $\epsilon$ is an unobservable random error with expected value $0$ and variance $\sigma^{2}\in(0,\infty)$. The random errors across different trials are assumed to be independent.

We remark that we restrict ourselves to linear homoscedastic regression models only for simplicity; the results of this paper can also be applied to nonlinear regression models through the approach of local optimality and to cases of non-constant error variance by means of a simple model transformation (e.g., [2], Chapters 17 and 22).

As is usual in the optimal design literature (e.g., Section 1.24 in [44]), we will formalize an approximate design as a probability measure $\xi$ on $\mathcal{X}$ with a finite support $\mathcal{S}(\xi):=\{\mathbf{x}_{1},\ldots,\mathbf{x}_{h}\}$, where $h$ is the support size, and $\mathbf{x}_{1},\ldots,\mathbf{x}_{h}\in\mathcal{X}$ are the support points. The interpretation is that for any $\mathbf{x}\in\mathcal{S}(\xi)$ the number $\xi(\mathbf{x}):=\xi(\{\mathbf{x}\})>0$ gives the proportion of the trials to be performed under the experimental conditions determined by the design point $\mathbf{x}$. For clarity, we will sometimes represent the design $\xi$ by the table

For a design $\xi$, the normalized information matrix is the positive semi-definite matrix

which is of size $m\times m$. We will assume that the problem at hand is nondegenerate in the sense that there is a design $\xi$ such that $\mathbf{M}(\xi)$ is nonsingular.

The quality of the design $\xi$ is measured by the “size” of $\mathbf{M}(\xi)$, assessed through a utility function known as a criterion. For simplicity,²²2It is straightforward to use the methods developed in this paper with any information function $\Phi$ that is strictly concave in the sense of Section 5.2 in [44]. we will focus on Kiefer’s criteria $\Phi_{p}$ with a parameter $p\in(-\infty,0]$ in their information-function versions, as defined in Chapter 6 of Pukelsheim [44]. For a singular $\mathbf{M}(\xi)$, the value of the $\Phi_{p}$-criterion is $0$, and for a non-singular $\mathbf{M}(\xi)$ it is defined by

This class includes the two most important criteria: $D$-optimality ($p=0$) and $A$-optimality ($p=-1$). We also implicitly cover the criterion of $I$-optimality (also called $V$- or $IV$-optimality), because the problem of $I$-optimal design is equivalent to the problem of $A$-optimal design in a linearly transformed model (e.g., Section 9.8 in Pukelsheim [44]).

For the chosen $p$, the $\Phi_{p}$-optimal design is any design $\xi^{*}$ that maximizes $\Phi_{p}(\mathbf{M}(\xi))$ in the set of all designs $\xi$. Note that our assumptions imply that there exists at least one $\Phi_{p}$-optimal design.

The most important result on optimal approximate designs is the so-called equivalence theorem; for Kiefer’s criteria it states that a design $\xi^{*}$ is $\Phi_{p}$-optimal if and only if $\mathbf{M}(\xi^{*})$ is non-singular and

with equality for all $\mathbf{x}\in\mathcal{S}(\xi^{*})$ (e.g., Section 7.20 in Pukelsheim [44]). In addition, for several standard linear regression models (e.g., those analyzed in Sections 3.2, 3.3 and 3.4), there exists a simple design $\xi^{*}$ such that $\mathbf{M}(\xi^{*})=\lambda\mathbf{I}_{m}$ for some $\lambda>0$ and

Then, by Corollary 4 of [26], the design $\xi^{*}$ is $\Phi_{p}$-optimal for all parameters $p$.

Let $p$ be fixed. While there can be an infinite number of $\Phi_{p}$-optimal designs, the properties of the Kiefer’s criteria and the assumption of nondegeneracy imply that the information matrix of all $\Phi_{p}$-optimal designs is the same; see, e.g., Section 8.14 in Pukelsheim [44]. We will denote this unique and necessarily nonsingular $\Phi_{p}$-optimal information matrix by $\mathbf{M}_{*}$. In the rest of this paper, if the criterion is general, fixed, or known from the context, we will drop the prefix $\Phi_{p}$-.

Of particular importance are the notions of a minimal optimal design and an absolutely minimal optimal design, introduced by Farrell et al. [16]. An optimal design $\xi^{*}$ is a minimal optimal design if there is no optimal design, whose support is a proper subset of $\mathcal{S}(\xi^{*})$.³³3Some authors use the expression “minimally-supported design” to denote any design with the support size equal to the number of parameters (Chang and Lin [7], Li and Majumdar [35] and many others), which is a different concept. An optimal design $\xi^{*}$ is called an absolutely minimal optimal design if there is no optimal design, whose support has fewer elements than $\mathcal{S}(\xi^{*})$.

As stated in the introduction, our aim is to characterize the set of all optimal designs, if we already know one optimal design. To this end, we will assume that we know a “maximal” optimal design $\xi^{*}$, which we define as an optimal design satisfying $\mathcal{S}(\xi^{*})\supseteq\mathcal{S}(\tilde{\xi}^{*})$ for any other optimal design $\tilde{\xi}^{*}$. We emphasize that in many models, such as those analysed in Section 3, the known and simplest optimal design $\xi^{*}$ satisfies

i.e, $\xi^{*}$ is maximal, owing to the equivalence theorem for $\Phi_{p}$-optimality.

Of course, for some optimal design problems, a maximal optimal design is not immediately apparent. Even then, it is usually possible to construct one by a procedure that inspects either the design points satisfying equality in (2), or the points of a theoretically derived finite set necessarily covering the supports of all optimal designs; see, e.g., Section 3 in [31] or Lemma 4.2 in [22].⁴⁴4Since there is a large variety of models with directly available maximal optimal designs to justify the broad applicability of our method, see Section 3, we decided not to detail such an inspection procedure or the theoretical results on possible supports of optimal designs in this paper. Note that some very specific models, such as the trigonometric regression on $[0,2\pi]$ (Section 9.16 in Pukelsheim [44]), do not admit any maximal optimal design. This is because in these models the union of the supports of all optimal designs is an infinite set and, according to the standard definition we use, any design only has a finite support. We exclude such optimal design problems from our further considerations.

While there can be many maximal optimal designs, their supports must be the same; we will denote the unique “maximum optimal support” by $\mathcal{Y}:=\{\mathbf{y}_{1},\ldots,\mathbf{y}_{d}\}$, where $d$ is the size of the maximum optimal support, and the regressor vectors corresponding to the trials at $\mathcal{Y}$ by $\mathbf{f}_{i}:=\mathbf{f}(\mathbf{y}_{i})$, for all $i=1,\ldots,d$.

Therefore, if our aim is to study only optimal designs, we can completely restrict our attention to the set $\mathcal{Y}$ and disregard the design points from $\mathcal{X}\setminus\mathcal{Y}$. Let us henceforth denote by $\Upsilon$ the set of all designs supported on $\mathcal{Y}$.

Because $\mathcal{Y}$ is finite, any design $\zeta$ in $\Upsilon$ can be characterized by the weights that it assigns to the design points $\mathbf{y}_{1},\ldots,\mathbf{y}_{d}$. We will call $\mathbf{w}_{\zeta}=(\zeta(\mathbf{y}_{1}),\ldots,\zeta(\mathbf{y}_{d}))^{\prime}$ the vector of weights of $\zeta$. This correspondence allows us to work with vectors from $\mathbb{R}^{d}$ instead of measures on $\mathcal{Y}$. The components of each $\mathbf{w}_{\zeta}$ must be non-negative and sum to one; the set of all weight-vectors on $\mathcal{Y}$ is therefore the probability simplex $\mathbb{P}:=\{\mathbf{w}\in\mathbb{R}^{d}:\mathbf{1}^{\prime}_{d}\mathbf{w}=1,\mathbf{w}\geq\mathbf{0}_{d}\}$. Analogously, given a vector of weights $\mathbf{w}\in\mathbb{P}$, the design $\zeta_{\mathbf{w}}$ determined by $\mathbf{w}$ satisfies $\zeta_{\mathbf{w}}(\mathbf{y}_{i})=w_{i}$ for all $i=1,\ldots,d$.

Because all optimal designs have the same information matrix, the set of all weight-vectors corresponding to the optimal designs is

where $\mathbf{M}(\zeta_{\mathbf{w}})=\sum_{i=1}^{d}w_{i}\mathbf{f}_{i}\mathbf{f}^{\prime}_{i}$. Since $\mathbb{P}$ is a bounded superset of $\mathbb{P}_{*}$, and the equalities $\mathbf{M}(\zeta_{\mathbf{w}})=\mathbf{M}_{*}$ are linear in $\mathbf{w}$, the set $\mathbb{P}_{*}$ is a polytope⁵⁵5Because the definition of the term “polytope” somewhat varies across the mathematical literature, we remark that by polytope we mean a bounded convex polyhedron, i.e., a bounded intersection of closed half-spaces. For an introduction to the theory of polytopes, see, e.g., Ziegler [55]. in $\mathbb{R}^{d}$. Therefore, we will say that $\mathbb{P}_{*}$ is the “polytope of optimal weights”. A compact expression of the polytope of optimal weights utilizes the matrix

where $\mathrm{vech}$ means the half-vectorization of a symmetric matrix. Therefore, the number of rows of A is $q:=\frac{1}{2}m(m+1)$.

The proof of Lemma 1 requires showing that $\mathbf{M}(\zeta_{\mathbf{w}})=\mathbf{M}_{*}$ for some $\mathbf{w}\geq 0$ entails $\mathbf{1}_{d}^{\prime}\mathbf{w}=1$; see Appendix B for the proofs. Note that (5) is geometrically the intersection of a finite number of halfspaces, that is, it is a H-representation of $\mathbb{P}_{*}$.

For the mathematical properties of the optimal design problem at hand, a key characteristic is the rank $s$ of $\mathbf{A}$, which we will call the “problem rank”. In the following lemma, “dimension” of a subset of a vector space refers to the affine dimension of its affine hull.

A simple consequence of Lemma 2 is that the polytope $\mathbb{P}_{*}$ is trivial ($t=0$), i.e., the optimal design is unique, if and only if $d=s$, which is if and only if $\mathbf{A}$ has full column rank. Note also that Lemma 2, parts (a) and (c), imply $t\leq d-m$. It follows that for the case $d=m$ the polytope $\mathbb{P}_{*}$ is trivial.

Because the dimension of $\mathbb{P}_{*}$ is $t<d$, $\mathbb{P}_{*}$ belongs to some proper affine subspace of $\mathbb{R}^{d}$. However, we can transform $\mathbb{P}_{*}$ to a full-dimensional polytope $\mathbb{T}^{*}\subset\mathbb{R}^{t}$ that has the same geometric properties (the number and incidence of the faces of all dimensions, the distance between vertices and so on), as follows.

Let $t>0$, let the columns of $\mathbf{U}=(\mathbf{u}_{1},\ldots,\mathbf{u}_{t})$ form an orthonormal basis of the null space $\mathcal{N}(\mathbf{A})$ and let $\xi^{*}$ be a maximal optimal design with the vector of weights $\mathbf{w}^{*}>\mathbf{0}_{d}$. Define an affine mapping $A$ at $\mathbb{R}^{t}$ by $A(\mathbf{t})=\mathbf{U}\mathbf{t}+\mathbf{w}^{*}$. Note that the image of $A$ is the affine space $\mathcal{A}=\mathcal{N}(\mathbf{A})+\mathbf{w}^{*}\supset\mathbb{P}_{*}$, and its inverse on $\mathcal{A}$ is $A^{-1}(\mathbf{w})=\mathbf{U}^{\prime}(\mathbf{w}-\mathbf{w}^{*})$. Therefore, $A$ is an affine bijection between $\mathbb{R}^{t}$ and $\mathcal{A}$, which means that the geometric structure of $\mathbb{P}_{*}$ is the same as the structure of

The interior of $\mathbb{T}_{*}$ contains the origin $\mathbf{0}_{t}$. The full-dimensional polytope $\mathbb{T}_{*}$ can be used for the numerical construction of specific optimal designs as outlined in Section 2.6, as well as for a visualization of the set of all optimal designs if $t\leq 3$.

Key elements of $\mathbb{P}_{*}$ are its vertices. Note that the number $\ell$ of these vertices is finite; we will denote them $\mathbf{v}_{1},\ldots,\mathbf{v}_{\ell}$. The representation (5) of $\mathbb{P}_{*}$ is in fact in the standard form of the feasible set of a linear program. As such, the well-established theory of linear programming can be used to give the following characterizations of $\mathbf{v}_{1},\ldots,\mathbf{v}_{\ell}$.

The central object of this paper is the set $\Upsilon_{*}$ of all $\Phi_{p}$-optimal designs, i.e.,

Here we formulate corollaries of the results of Section 2.3 in the standard optimal design terminology, so that it is clear how these results translate from the weights $\mathbf{w}\in\mathbb{R}^{d}$ to actual designs (i.e., finitely supported probability measures).

First, we show that the mathematical properties indeed translate to designs $\zeta\in\Upsilon$ by formalizing the transformation $\zeta\mapsto\mathbf{w}_{\zeta}$ as a linear isomorphism. Let $\Upsilon^{\pm}$ be the linear space of all signed measures supported on $\mathcal{Y}$. Clearly, $L:\zeta\to\mathbf{w}_{\zeta}=(\zeta(\mathbf{y}_{1}),\ldots,\zeta(\mathbf{y}_{d}))^{\prime}$ is a linear one-to-one mapping between $\Upsilon^{\pm}$ and $\mathbb{R}^{d}$, that is, a linear isomorphism, with inverse $L^{-1}:\mathbf{w}\to\zeta_{\mathbf{w}}$. Note that the image of $\Upsilon$ in $L$ is $\mathbb{P}$, and the image of $\Upsilon_{*}$ in $L$ is $\mathbb{P}_{*}$. The latter observation means that $\Upsilon_{*}$ is isomorphic to $\mathbb{P}_{*}$; we therefore call $\Upsilon_{*}$ the “polytope of optimal designs”.

The Carathéodory theorem (e.g., [47], Section 17) says that for any bounded set $\SS\subset\mathbb{R}^{r}$, every point in the convex hull of $\SS$ can be expressed as a convex combination of at most $r+1$ points of $\SS$. Taking Lemma 2 into account, we obtain:

Theorem 2 is a strengthening of the standard “Carathéodory theorem” used in optimal design (e.g., Theorem 2.2.3 in Fedorov [17] or Section 8.2. in Pukelsheim [44]): we obtain the standard result for $s=q$, but in general $s\leq q$.

Recall that $\zeta_{\mathbf{w}}\in\Upsilon$ is the design determined by the weights $\mathbf{w}\in\mathbb{P}$ and $\mathbf{v}_{1},\ldots,\mathbf{v}_{\ell}$ are the vertices of $\mathbb{P}_{*}$. We will call the corresponding designs $\nu_{j}=\zeta_{\mathbf{v}_{j}}\in\Upsilon_{*}$, $j=1,\ldots,\ell$, the “vertex optimal designs” (VODs). In simpler terms, for each vertex $\mathbf{v}_{j}=(v_{j,1},\ldots,v_{j,d})^{\prime}$ of $\mathbb{P}_{*}$, $j\in\{1,\ldots,\ell\}$, the uniquely determined vertex optimal design is

Since the mapping $\mathbf{w}\mapsto\zeta_{\mathbf{w}}$ is a linear isomorphism between $\mathbb{P}_{*}$ and $\Upsilon_{*}$, the role of the vertices $\mathbf{v}_{1},\ldots,\mathbf{v}_{\ell}$ in $\mathbb{P}_{*}$ is analogous to the role of the VODs $\nu_{1},\ldots,\nu_{\ell}$ in $\Upsilon_{*}$.

As a consequence of Lemma 2 we see that the dimension of $\Upsilon_{*}$, understood as the dimension of its affine hull in $\Upsilon^{\pm}$, is $t=d-s$. The knowledge of $t$ and the inequality $s\geq m$ provides another variant of the Carathéodory theorem expressed in terms of the VODs:

The main theorem characterizing the VODs, based directly on Theorem 1, is

Note that Theorem 4, specifically part $1\Leftrightarrow 3$, implies that the support size of any VOD is between $m$ and $s$. Theorem 4 also provides bounds on the number $\ell$ of the VODs in terms of the numbers $d$, $s$ and $t$; see Appendix D.

In Examples 3 to 5, we provided the complete characterization of the optimal design polytope via its vertices for Model (1). However, such analytical derivations become tedious with increasing number of factors and more complex models.

Fortunately, the problems of finding all VODs for many standard models (including Model (1)) are “rational”, that is, the polytope of optimal weights has a rational H-representation (5). More precisely, it means that the matrices $\mathbf{M}_{*}$ and $\mathbf{f}_{i}\mathbf{f}^{\prime}_{i}$ ($i=1,\ldots,d$) have rational elements. These conditions are automatically satisfied if $\mathbf{f}_{1},\ldots,\mathbf{f}_{d}$ have rational elements, and the weights of at least one optimal design are rational numbers.

For a rational optimal design problem, all VODs $\nu_{1},\ldots,\nu_{\ell}$ have rational weights. Moreover, a rational optimal design problem allows for the use of error-free rational-arithmetic algorithms to construct the theoretical forms of all VODs. In Section 3, we therefore employ computer-assisted construction of VODs that are (analytically) optimal, rather than the manual construction of such designs.

The rational form of the weights of all VODs means that they are exact designs of some sizes $N_{1},\ldots,N_{\ell}$. It is of theoretical interest that in the rational optimal design problems, any optimal approximate design is a convex combination of at most $t+1$ exact optimal designs; see Theorem 3. From the more practical point of view, the VODs can be combined to provide an entire pool of optimal designs that can be realized with a finite number of trials. In particular, for some $N_{0}$ we can achieve an exact design of any size $N\geq N_{0}$ that is a multiple of the greatest common divisor $\textrm{gcd}(N_{1},\ldots,N_{\ell})$, cf. Nijenhuis and Wilf [40]. We illustrate this idea for all models analyzed in Section 3.

Finally, note that for a rational design problem, its rank $s$ can also be calculated in a perfectly reliable way using the rational arithmetic, even without the need to calculate the VODs. The problem rank alone provides the dimension of the set of all information matrices on $\mathcal{Y}$ (Lemma 2), an upper bound on the minimal support size (Theorem 2) and the dimension $t=d-s$ of the polytope $\Upsilon_{*}$ of optimal designs (Lemma 2).

Suppose that the optimal design is not unique, and we have the complete list of all VODs, characterizing the set $\Upsilon_{*}$ of all optimal designs. Then we can select one or more optimal designs with special properties, to broaden the choice for the experimenter.

Consider the model

The case of a single factor ($k=1$) is trivial and for two factors ($k=2$), all $\Phi_{p}$-optimal designs are unique; see Section 8.6 in Pukelsheim [44].

Let $k\geq 3$. For $j=1,\ldots,k$, let $\mathcal{V}(j)$ be the set of all vectors in $\{0,1\}^{k}$ with exactly $j$ components equal to $1$. The following claims are straightforward to verify via (2) and (4); cf. Pukelsheim [44], page 376: If $k$ is odd, the maximum optimal support is $\mathcal{Y}_{D}=\mathcal{Y}_{A}=\mathcal{V}(k/2+1/2)$ for both $D$- and $A$-optimality. If $k$ is even, the maximum optimal support is $\mathcal{Y}_{D}=\mathcal{V}(k/2)\cup\mathcal{V}(k/2+1)$ for $D$-optimality and $\mathcal{Y}_{A}=\mathcal{V}(k/2)$ for $A$-optimality. For any $k$, the design uniform on $\mathcal{Y}_{D}$ is maximal $D$-optimal, and the design uniform on $\mathcal{Y}_{A}$ is maximal $A$-optimal.

Consequently, it is possible to apply rational arithmetic calculations to prove that for $k\leq 5$ the uniform optimal design on $\mathcal{Y}_{D}$ is the unique $D$-optimal design, and the uniform design on $\mathcal{Y}_{A}$ is the unique $A$-optimal design.

For $k=6$ and the $D$-criterion, the dimension of $\Upsilon_{*}$ is $t=14$, and there are $\ell=150$ VODs. To describe the set of all VODs, consider two designs to be isomorphic, if they are the same up to a relabeling of the factors. Then, the VODs can be split into two classes of isomorphic designs. The first class consists of the $30$ VODs that are isomorphic to

All $30$ VODs in this class can be obtained by permuting the (first $6$) rows of the table. Note that permuting the columns corresponds to the same VOD. One of these VODs has been computed by integer quadratic programming as an exact $D$-optimal design of size $7K$, $K\in\mathbb{N}$; see Table 1 in Filová and Harman [18]. For $k=6$, these $30$ VODs are the absolutely minimal $D$-optimal designs.

The second class of VODs for $D$-optimality and $k=6$ consists of the $120$ VODs isomorphic to the uniform design

For $k=6$ and for the criterion of $A$-optimality, $\Upsilon_{*}$ has dimension $t=5$, and there are $\ell=12$ VODs, all isomorphic to

One of these VODs has been identified by integer quadratic programming as an $A$-optimal exact design of size $10K$, $K\in\mathbb{N}$; see Table 2 in Filová and Harman [18]. For $k=6$, these $12$ VODs are the absolutely minimal $A$-optimal designs. A $3$-dimensional projection of $\mathbb{P}_{*}$ and its $12$ vertices are depicted in Figure 1.

For $k=7$, the sets of $D$- and $A$-optimal designs coincide. The dimension of $\Upsilon_{*}$ is $t=14$ and it has $\ell=150$ vertices; their structure is similar to the case of $D$-optimality for $k=6$.

For $k=8$ and the criterion of $D$ optimality, it holds that the size of the maximum optimal support is $d=126$, the problem rank is $s=36$, and the dimension of $\Upsilon_{*}$ is $t=90$. For $k=8$ and the criterion of $A$-optimality we have $d=70$, $s=38$ and $t=42$. We were not able to enumerate the individual VODs for this size.

Relations to block designs. For Model (10), a $k\times b$ matrix of support points of a design on $\mathcal{Y}_{D}$ or $\mathcal{Y}_{A}$ can be viewed as the incidence matrix of a binary block design with $k$ treatments and $b$ blocks. In this view, the support points of the $30$ VODs isomorphic to (11) are partially balanced designs (PBDs; see Wilson [54]) with $r=4$, $\lambda=2$, and $b=7$. Here, $r$ is the common replication number, and $\lambda$ the common concurrence. Similarly, the $120$ VODs isomorphic to (12) are PBDs with $r=12$, $\lambda=6$, and $b=21$. In addition, the support points of (13) correspond to a balanced incomplete block design (which is a PBD with a constant block size) for $6$ treatments, $r=5$, $\lambda=2$, $b=10$, and a common block size $3$ (cf. Pukelsheim [44], p. 378, Cheng [8]).

In this section we consider the model