Improvements on removing non-optimal support points in D-optimum design algorithms

Introduction

denote the information matrix. Suppose that there exists a design with nonsingular information matrix and let $\Xi^{+}$ be the set of such designs. Let $\xi^{*}$ denote a $D$-optimum design, that is, a measure in $\Xi$ that maximizes $\det\mathbf{M}(\xi)$, see, e.g., (Fedorov 1972). Note that a $D$-optimum design always exists and that the $D$-optimum information matrix $\mathbf{M}_{*}=\mathbf{M}(\xi^{*})$ is unique. For any $\xi\in\Xi^{+}$ denote $d(\xi,\cdot):\mathcal{X}\to[0,\infty)$ the variance function defined by

The celebrated Kiefer-Wolfowitz Equivalence Theorem (1960) writes as follows.

The following three statements are equivalent:

$\max_{\mathbf{x}\in{\mathcal{X}}}d(\xi^{*},\mathbf{x})=m$;

$\xi^{*}$ minimizes $\max_{\mathbf{x}\in{\mathcal{X}}}d(\xi,\mathbf{x})$, $\xi\in\Xi^{+}$.

Hence, (ii) of Theorem 1 implies that for any support point $\mathbf{x}_{*}$ of the design $\xi^{*}$ (i.e., for a point satisfying $\xi^{*}(\mathbf{x}_{*})>0$), we have

In the next section we show that the equality (1) can be used to prove that

where $\lambda_{1}^{*}$ depends on $\xi$ only via the maximum of $d(\xi,\cdot)$ over the design space ${\mathcal{X}}$. Hence, we can test candidate support points by using any finite number of design measures $\xi\in\Xi^{+}$, e.g., those that are generated by a design algorithm on its way towards the optimum: any point that does not pass the test defined by $\xi^{k}$ of iteration $k$ need not be considered for further investigations and can thus be removed from the design space.

A necessary condition for candidate support points

For $\xi$ a design in $\Xi^{+}$ denote $\mathbf{M}=\mathbf{M}(\xi)$,

and $\lambda_{1}\leq\lambda_{2}\leq\cdots\leq\lambda_{m}$ the eigenvalues of $\mathbf{H}$. Notice that $\lambda_{1}>0$ and that the eigenvalues depend on the design $\xi$ as well as on the $D$-optimum information matrix $\mathbf{M}_{*}$. Let $\mathbf{x}_{*}$ be a support point of a $D$-optimum design and let $\mathbf{y}_{*}=\mathbf{H}^{-1/2}\mathbf{M}^{-1/2}\mathbf{x}_{*}$. The equality (1) can be written in the form $\mathbf{y}_{*}^{\top}\mathbf{y}_{*}=m$ which implies:

To be able to use the inequality (2), we need to derive a lower bound $\lambda_{1}^{*}$ on $\lambda_{1}$ that does not depend on the unknown matrix $\mathbf{M}_{*}$.

For $m=1$ we directly obtain the lower bound $\lambda_{1}\geq\lambda_{1}^{*}=1$. For $m>1$, the Lagrangian for the minimisation of $\lambda_{1}$ subject to $\sum_{i=1}^{m}\lambda_{i}^{-1}\leq m$ and $\sum_{i=1}^{m}\lambda_{i}\leq m+\epsilon$ is given by

with $\lambda=(\lambda_{1},...,\lambda_{m})^{\top}$, $\mu_{1},\mu_{2}\geq 0$. The stationarity of $\mathcal{L}(\lambda,\mu_{1},\mu_{2})$ with respect to the $\lambda_{i}$’s and the Kuhn-Tucker conditions

give $\lambda_{i}=L$ for $i=2,...,m$, with $\lambda_{1}$ and $L$ satisfying

and $\lambda_{i}^{*}=L^{*}=(m-1)/(m-1/\lambda_{1}^{*})\geq 1$, $i=2,...,m$. Notice that the bound (4) gives $\lambda_{1}^{*}=1$ when $m=1$ and can thus be used for any dimension $m\geq 1$. By substituting $\lambda_{1}^{*}$ for $\lambda_{1}$ in (2) we obtain the following result.

For any design $\xi\in\Xi^{+}$, any point $\mathbf{x}_{*}\in\mathcal{X}$ such that

where $\epsilon=\max_{\mathbf{x}\in\mathcal{X}}d(\xi,\mathbf{x})-m$, cannot be a support point of a D-optimum design measure.

The uniform probability measure $\eta$ on $\mathbf{y}_{1},...,\mathbf{y}_{k}$ is $D$-optimum on $\mathcal{X}/\{\mathbf{x}_{*}\}$, as can be directly verified by checking (ii) of the Equivalence Theorem 1. On the other hand, $\eta$ is not $D$-optimum on $\mathcal{X}$ since $\mathbf{x}_{*}^{\top}\mathbf{M}^{-1}(\eta)\mathbf{x}_{*}=b\,m>m$, which implies that $\mathbf{x}_{*}$ must support a $D$-optimum design on $\mathcal{X}$.

Figure 1 presents a typical evolution of $q(k)$ as a function of $\log(k)$ for $\xi^{0}$ uniform on ${\mathcal{X}}$ and shows the superiority of the test (5) over (6). The improvement is especially important in the first iterations, when the design $\xi^{k}$ is far from the optimum. Define $k^{*}(\delta)$ as the number of iterations required to reach a given precision $\delta$,

with $\epsilon(\xi^{k})$ defined by (3). Notice that from the concavity of $\log\det\mathbf{M}(\xi)$ we have

Table 1 shows the influence on the algorithm (7) of the cancellation of points based on the tests (5) and (6), in terms of $k^{*}(\delta)$, of the corresponding computing time $T(\delta)$, the number of support points $n(\delta)$ of $\xi^{k^{*}(\delta)}$ and the first iteration $k_{10}$ when $\xi^{k}$ has 10 support points or less, with $\delta=10^{-3}$. The results are averaged over 1000 independent problems. The values of $k^{*}(\delta)$ and $k_{10}$ are rounded to the nearest larger integer, the computing time for the algorithm with the cancellation of points based on (5) is taken as reference and set to 1 (the algorithm without cancellation was at least 4.5 times slower in all the 1000 repetitions). Although cancelling points has little influence on the number of iterations $k^{*}(\delta)$, is renders the iterations simpler: on average the introduction of the test (5) in the algorithm (7) makes it about 30 times faster.

The influence of the cancellation on the performance of the algorithm can be further improved as follows. Let $(k_{j})_{j}$ denote the subsequence corresponding to the iterations where some points are removed from ${\mathcal{X}}$. We have $j\leq q(0)$, the cardinality of the initial ${\mathcal{X}}$, and the convergence of the algorithm (7) is therefore maintained whatever the heuristic rule used at the iterations $k_{j}$ for updating the weights of the points that stay in ${\mathcal{X}}$ (provided these weights remain strictly positive). The following one has been found particularly efficient on a series of examples: for all $t\in T_{j}$, the set of indices corresponding to the points that stay in ${\mathcal{X}}$ at iteration $k_{j}$, replace $w_{t}^{k_{j}}$ by

for some $A\geq 1$. A final remark is that by including the test (5) in the algorithm (7) one can in general quickly identify potential support points for an optimum design. When the number $n$ of these points is small enough, switching to a more standard convex-programming algorithm for the optimization of the $n$ associated weights might then form a very efficient strategy.

References