Monotonic convergence of a general algorithm for computing optimal designs

Yaming Yu

A general class of algorithms

Optimal experimental design (approximate theory) is a well-developed area, and we refer to Kiefer (Kiefer74 1974), Silvey (Silvey80 1980), Pázman (Paz86 1986) and Pukelsheim (Pu93 1993) for a general introduction and basic results. We consider computational aspects of optimal designs, focusing on a finite design space $\mathcal{X}=\{x_{1},\ldots,x_{n}\}$ . Suppose the probability density or mass function of the response is specified as $p(y|x,\theta)$ where $\theta=(\theta_{1},\ldots,\theta_{m})^{\top}$ is the parameter of interest. Let $A_{i}$ denote the $m\times m$ expected Fisher information matrix from a unit assigned to $x_{i}$ with the $(j,k)$ entry [the expectation is with respect to $p(y|x_{i},\theta)$ ]

The moment matrix, as a function of the design measure $w=(w_{1},\ldots,w_{n})$ , is defined as

which is proportional to the Fisher information for $\theta$ when the number of units assigned to $x_{i}$ is proportional to $w_{i}$ . Here $w\in\bar{\Omega}$ , and $\bar{\Omega}$ denotes the closure of $\Omega=\{w\colon\;w_{i}>0,\sum_{i=1}^{n}w_{i}=1\}$ . Throughout we assume that $A_{i}$ are well defined and hence nonnegative definite. The set

is assumed nonempty. Our approach may conceivably extend to the case where $M(w)$ is allowed to be singular, by using generalized inverses, although we do not pursue this here.

Given an optimality criterion $\phi$ defined on positive definite matrices, the goal is to maximize $\phi(M(w))$ with respect to $w\in\Omega_{+}$ . Typical optimality criteria include:

the D-criterion $\phi_{0}(M)=\log\det(M)$ ,

the A-criterion $\phi_{-1}(M)=-\operatorname{tr}(M^{-1})$ ,

more generally, the $p$ th mean criterion $\phi_{p}(M)=-\operatorname{tr}(M^{p}),p<0$ and

the c-criterion $\phi_{-1,c}(M)=-c^{\top}M^{-1}c$ , where $c$ is a nonzero constant vector.

Often only a linear combination $K^{\top}\theta$ , for example, a subvector of $\theta$ , is of interest. The Fisher information for $K^{\top}\theta$ is naturally defined as $(K^{\top}M^{-1}K)^{-1}$ , assuming invertibility [Pukelsheim (Pu93 1993)]. We may therefore consider the D- and A-criteria for $K^{\top}\theta$ defined, respectively, as

The c-criterion is a special case of $\phi_{-1,K}(M)$ . Motivations for such optimality criteria are well known. In a linear problem, the A-criterion seeks to minimize the sum of variances of the best linear unbiased estimators (BLUEs) for all coordinates of $\theta$ while the c-criterion seeks to minimize the variance of the BLUE for $c^{\top}\theta$ . Similar interpretations (with asymptotic arguments) apply to nonlinear problems.

In general $M(w)$ also depends on the unknown parameter $\theta$ which complicates the definition of an optimality criterion. A simple solution is to maximize $\phi(M(w))$ with $\theta$ fixed at a prior guess $\theta^{*}$ ; this leads to local optimality [Chernoff (Chernoff53 1953)]. Local optimality may be criticized for ignoring uncertainty in $\theta$ . However, in a situation where real prior information is available, or where the dependence of $M$ on $\theta$ is weak, it is nevertheless a viable approach and has been adopted routinely [see, e.g., Li and Majumdar (LM08 2008)]. Henceforth we assume a fixed $\theta^{*}$ and suppress the dependence of $M$ on $\theta$ . Possible extensions are mentioned in Section 5.

Optimal designs do not usually come in closed form. As early as Wynn (Wynn72 1972), Fedorov (Fed72 1972), Atwood (Atwood73 1973) and Wu and Wynn (WW78 1978), and as late as Torsney (Torsney07 2007), Harman and Pronzato (HP07 2007) and Dette, Pepelyshev and Zhigljavsky (DPZ08 2008), various procedures have been studied for numerical computation. We shall focus on the following multiplicative algorithm [Titterington (Titterington76 1976, Titterington78 1978), Silvey, Titterington and Torsney (STT78 1978)] which is specified through a power parameter $\lambda\in(0,1]$ .

Set $\lambda\in(0,1]$ and $w^{(0)}\in\Omega$ . For $t=0,1,\ldots,$ compute

For a heuristic explanation, observe that (2) is equivalent to

The value of $\partial\phi(M(w))/\partial w_{i}$ indicates the amount of gain in information, as measured by $\phi$ , by a slight increase in $w_{i}$ , the weight on the $i$ th design point. So (3) can be seen as adjusting $w$ so that relatively more weight is placed on design points whose increased weight may result in a larger gain in $\phi$ . If $\phi$ is increasing and concave, then a convenient convergence criterion, based on the general equivalence theorem [Kiefer and Wolfowitz (KW60 1960), Whittle (Whittle73 1973)], is

where $\bar{d}(w)\equiv\sum_{i=1}^{n}w_{i}d_{i}(w)$ and $\delta$ is a small positive constant.

Algorithm I is remarkable in its generality. For example, little restriction is placed on the underlying model $p(y|x,\theta)$ . Part of the reason, of course, is that we focus on Fisher information and local optimality, which essentially reduces the problem to a linear one.

There exists a large literature on Algorithm I and its relatives [see, e.g., Titterington (Titterington76 1976, Titterington78 1978), Silvey, Titterington and Torsney (STT78 1978), Pázman (Paz86 1986), Fellman (Fellman89 1989), Pukelsheim and Torsney (PT91 1991), Torsney and Mandal (TM06 2006), Harman and Pronzato (HP07 2007), Dette, Pepelyshev and Zhigljavsky (DPZ08 2008) and Torsney and Martín-Martín (TM09 2009)]. One feature that has attracted much attention is that Algorithm I appears to be monotonic, that is, $\phi(M(w^{(t)}))$ increases in $t$ , at least in some special cases. For example, when $\phi=\phi_{0}$ (for D-optimality) and $\lambda=1$ , Titterington (Titterington76 1976) and Pázman (Paz86 1986) have shown monotonicity using clever probabilistic and analytic inequalities [see also Dette, Pepelyshev and Zhigljavsky (DPZ08 2008) and Harman and Trnovská (HT09 2009)]. Algorithm I is also known to be monotonic for $\phi=\phi_{-1,K}$ as in (1), assuming $\lambda=1/2$ and $A_{i}$ are rank-one [Fellman (Fellman74 1974), Torsney (Torsney83 1983)]. Monotonicity is important because convergence then holds under mild assumptions (see Section 4). Results in these special cases suggest a monotonic convergence theory for a broad class of $\phi$ which is also supported by numerical evidence presented in some of the references above.

Main result

We aim to state general conditions on $\phi$ that ensure that Algorithm I converges monotonically. As a consequence certain known theoretical results are unified and generalized, and one particular conjecture [Titterington (Titterington78 1978)] is confirmed. Define

The functions $\phi$ and $\psi$ are assumed to be differentiable on invertible matrices. Our conditions are conveniently stated in terms of $\psi$ . As usual, for two symmetric matrices, $M_{1}\leq(<)M_{2}$ means $M_{2}-M_{1}$ is nonnegative (positive) definite.

or, equivalently, $\psi^{\prime}(M)$ is nonnegative definite for positive definite $M$ .

for $\alpha\in,M_{1},M_{2}>0$ . Equivalently,

Condition (5) is usually satisfied by any reasonable information criterion[Pukelsheim (Pu93 1993)]. Also note that, if (5) fails, then $\partial\phi(M(w))/\partial w_{i}$ on the right-hand side of (3) is not even guaranteed to be nonnegative. The real restriction is the concavity condition (6). For example, (6) is not satisfied by $\psi_{p}(M)=-\phi_{p}(M^{-1})$ (the $p$ th mean criterion) when $p<-1$ . [It is usually assumed that $\phi(M)$ , rather than $\psi(M)$ , is concave.] Nevertheless, (6) is satisfied by a wide range of criteria, including the commonly used D-, A- or c-criteria [see cases (i) and (ii) in the illustration of the main result below].

Assume (5) and (6). Assume that in iteration (2), with $0<\lambda\leq 1$ , we have

In other words, under mild conditions which ensure that (2) is well defined [specifically, the denominator in (2) is nonzero], (5) and (6) imply that (2) never decreases the criterion $\phi$ . Let us illustrate Theorem 1 with some examples. For simplicity, in (i)–(iv) we display formulae for $\lambda=1$ only, although monotonicity holds for all $\lambda\in(0,1]$ .

Then $\psi_{p}(M)\equiv-\phi_{p}(M^{-1})$ satisfies (5) and (6). By Theorem 1, Algorithm I is monotonic for $\phi=\phi_{p},p\in$ . This generalizes the previously known cases $p=0$ and $p=-1$ (with particular values of $\lambda$ ). The iteration (2) reads

More generally, given a full rank $m\times r$ matrix $K$ ( $r\leq m$ ), consider

Then $\psi_{p,K}(M)$ satisfies (5) and (6). By Theorem 1, Algorithm I is monotonic for $\phi=\phi_{p,K},p\in$ . The iteration (2) reads

In particular, taking $r=1,K=c$ (an $m\times 1$ vector) and $p=-1$ in case (ii), we obtain that Algorithm I is monotonic for the c-criterion $\phi_{-1,c}$ . The iteration (8) reduces to

As noted by a referee, with $p=-1$ , the choice $\lambda=1$ may lead to an oscillating behavior in the sense that $w^{(t)}$ alternates between two points at which $\phi_{-1,c}(M(w))$ takes the same value. While this does not contradict Theorem 1, it suggests that other values of $\lambda$ are more desirable for fast convergence. Following Fellman (Fellman74 1974) and Torsney (Torsney83 1983), a practical recommendation is $\lambda=1/2$ in the $p=-1$ case.

Consider another example of case (ii), with $p=0,r=m-1$ and $K=(0_{r},I_{r})^{\top}$ . Henceforth $0_{r}$ denotes the $r\times 1$ vector of zeros, and $I_{r}$ denotes the $r\times r$ identity matrix. Assume $A_{i}=x_{i}x_{i}^{\top},x_{i}^{\top}=(1,z_{i}^{\top})$ and $z_{i}$ is $(m-1)\times 1$ . This corresponds to a D-optimal design problem for $(\theta_{2},\ldots,\theta_{m})$ under the linear model,

where the parameter is $\theta=(\theta_{1},\theta_{2},\ldots,\theta_{m})^{\top}$ . That is, interest centers on all coefficients other than the intercept. Nevertheless, as far as the design measure $w$ is concerned, the optimality criterion, $\phi_{0,K}(M)$ , coincides with $\phi_{0}(M)$ , that is,

Thus (9) satisfies $\det M(w^{(t+1)})\geq\det M(w^{(t)})$ .

Monotonicity of (9) has been conjectured since Titterington (Titterington78 1978), and considerable numerical evidence has accumulated over the years. Recently, extending the arguments of Pázman (Paz86 1986), Dette, Pepelyshev and Zhigljavsky (DPZ08 2008) have obtained results which come very close to resolving Titterington’s conjecture. Nevertheless, we have been unable to extend their arguments further. Instead we prove the general Theorem 1 using a different approach, and settle this conjecture as a consequence.

The proof of Theorem 1 is achieved by using a method of auxiliary variables. When a function $f(w)$ [e.g., $-\det M(w)$ ] to be minimized is complicated, we introduce a new variable $Q$ and a function $g(w,Q)$ such that $\min_{Q}g(w,Q)=f(w)$ for all $w$ , thus transforming the problem into minimizing $g(w,Q)$ over $w$ and $Q$ jointly. Then we may use an iterative conditional minimization strategy on $g(w,Q)$ . This is inspired by the EM algorithm [Dempster, Laird and Rubin (DLR77 1977), Meng and van Dyk (MV97 1997); in particular, see Csiszár and Tusnady’s (CT84 1984) interpretation; see Yu (Yu08 2008) for a related interpretation of the data augmentation algorithm].

In Section 3 we analyze Algorithm I using this strategy. Although attention is paid to the mathematics, our focus is on intuitively appealing interpretations which may lead to further extensions of Algorithm I with the same desirable monotonicity properties. If the algorithm is monotonic, then convergence can be established under mild conditions (Section 4). Section 5 contains an illustration with optimal designs for a simple logistic regression model.

Explaining the monotonicity

A key observation is that the problem of maximizing $\phi(M(w))$ , or, equivalently, minimizing $\psi(M^{-1}(w))$ can be formulated as a joint minimization over both the design and the estimator. Specifically, let us compare the original Problem P1 with its companion Problem P2. Throughout $A^{1/2}$ denotes the symmetric nonnegative definite (SNND) square root of an SNND matrix $A$ .

Minimize $-\phi(M(w))\equiv\psi((\sum_{i=1}^{n}w_{i}A_{i})^{-1})$ over $w\in\Omega$ .

over $w\in\Omega$ and $Q$ [an $m\times(mn)$ matrix], subject to $QG=I_{m}$ , where

Though not immediately obvious, Problems P1 and P2 are equivalent, and this may be explained in statistical terms as follows. In (10), $Q\Delta_{w}Q^{\top}$ is simply the variance matrix of a linear unbiased estimator, $QY$ , of the $m\times 1$ parameter $\theta$ in the model

where $Y$ is the $(mn)\times 1$ vector of observations. The constraint $QG=I_{m}$ ensures unbiasedness. [Note that $G$ is full-rank since $M(w)$ is nonsingular by assumption.] Of course, the weighted least squares (WLS) estimator is the best linear unbiased estimator, having the smallest variance matrix (in the sense of positive definite ordering) and, by (5), the smallest $\psi$ for that matrix. It follows that, for fixed $w$ , $g(w,Q)$ is minimized by choosing $QY$ as the WLS estimator,

That is, Problem P2 reduces to Problem P1 upon minimizing over $Q$ .

Since Problem P2 is not immediately solvable, it is natural to consider the subproblems: (i) minimizing $g(w,Q)$ over $Q$ for fixed $w$ and (ii) minimizing $g(w,Q)$ over $w$ for fixed $Q$ . Part (ii) is again formulated as a joint minimization problem. For a fixed $m\times(mn)$ matrix $Q$ such that $QG=I_{m}$ , let us consider Problems P3 and P4.

Minimize $g(w,Q)$ as in (10) over $w\in\Omega$ .

over $w\in\Omega$ and the $m\times m$ positive-definite matrix $\Sigma$ .

The concavity assumption (7) implies that

with equality when $\Sigma=Q\Delta_{w}Q^{\top}$ , that is, Problem P4 reduces to Problem P3 upon minimizing over $\Sigma$ .

Since Problem P4 is not immediately solvable, it is natural to consider the subproblems: (i) minimizing $h(\Sigma,w,Q)$ over $\Sigma$ for fixed $w$ and $Q$ and (ii) minimizing $h(\Sigma,w,Q)$ over $w$ for fixed $\Sigma$ and $Q$ . Part (ii), which amounts to minimizing

admits a closed-form solution: if we write $Q=(Q_{1},\ldots,Q_{n})$ where each $Q_{i}$ is $m\times m$ , then $w_{i}^{2}$ should be proportional to $\operatorname{tr}(Q_{i}^{\top}\psi^{\prime}(\Sigma)Q_{i})$ . But Algorithm I may not perform an exact minimization here [see (16)].

Based on the above discussion, we can express Algorithm I as an iterative conditional minimization algorithm involving $w,Q$ and $\Sigma$ . At iteration $t$ , define

The choice of $w^{(t+1)}$ leads to (16) as follows. After simple algebra, the iteration (2) becomes

Since $0<\lambda\leq 1$ , Jensen’s inequality yields

which produces (16). Choosing $\lambda=1/2$ , that is, $w_{i}^{(t+1)}\propto\sqrt{r_{i}}$ , leads to exact minimization in (16); choosing $\lambda=1$ yields equality in (16). But any choice of $w^{(t+1)}$ that decreases $h(\Sigma^{(t)},w,Q^{(t)})$ at (16) would have resulted in the desired inequality,

We may allow $\lambda$ to change from iteration to iteration, and monotonicity still holds, as long as $\lambda\in(0,1]$ . See Silvey, Titterington and Torsney (STT78 1978) and Fellman (Fellman89 1989) for investigations concerning the choice of $\lambda$ . Also note that we assume $w_{i}^{(t)},w_{i}^{(t+1)}>0$ for all $i$ . This is not essential, however, because (i) the possibility of $w_{i}^{(t)}=0$ can be handled by restricting our analysis to all design points $i$ such that $w_{i}^{(t)}>0$ , and (ii) the possibility of $w_{i}^{(t+1)}=0$ can be handled by a standard limiting argument. Monotonicity holds as long as $M(w^{(t)})$ and $M(w^{(t+1)})$ are both positive definite, as noted in the statement of Theorem 1.

Global convergence

Monotonicity (Theorem 1) plays an important role in the following convergence theorem.

(b) Assume (2) is strictly monotonic, that is,

(d) Assume that, if $M$ (a positive definite matrix) tends to $M^{*}$ such that $\phi(M)$ increases monotonically, then $M^{*}$ is nonsingular.

Let $w^{(t)}$ be generated by (2) with $w^{(0)}_{i}>0$ for all $i$ . Then:

all limit points of $w^{(t)}$ are global maxima of $\phi(M(w))$ on $\Omega_{+}$ , and

as $t\to\infty$ , $\phi(M(w^{(t)}))$ increases monotonically to $\sup_{w\in\Omega_{+}}\phi(M(w))$ .

The proof of Theorem 2 is somewhat subtle. Standard arguments show that all limit points of $w^{(t)}$ are fixed points of the mapping $T$ . This alone does not imply convergence to a global maximum, however, because there often exist sub-optimal fixed points on the boundary of $\Omega$ . (Global maxima occur routinely on the boundary also.) Our goal is therefore to rule out possible convergence to such sub-optimal points; details of the proof are presented in Yu (Yu09 2009), an extended version of this paper. We shall comment on conditions (a)–(d).

Condition (a) ensures that starting with $w^{(0)}\in\Omega_{+}$ , all iterations are well defined. Moreover, if $w^{(0)}_{i}>0$ for all $i$ , then $w^{(t)}_{i}>0$ for all $t$ and $i$ . This highlights the basic idea that, in order to converge to a global maximum $w^{*}$ , the starting value $w^{(0)}$ must assign positive weight to every support point of $w^{*}$ . Such a requirement is not necessary for monotonicity. On the other hand, assigning weight to nonsupporting points of $w^{*}$ tends to slow the algorithm down. Hence methods that quickly eliminate nonoptimal support points are valuable [Harman and Pronzato (HP07 2007)].

Condition (b) simply says that unless $w$ is a fixed point, the mapping $T$ should produce a better solution. Let us assume (5), (7) and condition (a) so that Theorem 1 applies. Then, by checking the equality condition in (16), it is easy to see that condition (b) is satisfied if $0<\lambda<1$ . [The argument leading to (19) technically assumes that all coordinates of $w$ are nonzero, but we can apply it to the appropriate subvector of $w$ .] If $\lambda=1$ , then (16) reduces to an equality. However, by checking the equality conditions in (17) and (18), we can show that condition (b) is satisfied if $\psi$ is strictly increasing and strictly concave:

Conditions (c) and (d) are technical requirements that concern $\phi$ alone. Condition (c) ensures uniqueness of the optimal moment matrix which simplifies the analysis. Condition (d) ensures that positive definiteness of $M(w)$ is maintained in the limit. Conditions (c) and (d) are satisfied by $\phi=\phi_{p}$ with $p\leq 0$ , for example.

Let us mention a typical example of Theorem 2.

Assume $A_{i}\neq 0,w^{(0)}_{i}>0,i=1,\ldots,n$ , and $M(w^{(0)})>0$ . Then the conclusion of Theorem 2 holds for Algorithm I with $\phi=\phi_{0}$ .

Conditions (a), (c) and (d) are readily verified. Condition (b) is satisfied by (4) and (4). The claim follows from Theorem 2.

When (4) or (4) fails, and $\lambda=1$ , it is often difficult to appeal to Theorem 2 because strict monotonicity [condition (b)] may not hold. We illustrate this with an example where the monotonicity is not strict, and the algorithm does not converge [see Pronzato, Wynn and Zhigljavsky (PrWyZh2000 2000), Chapter 7; also the remark in case (iii) following Theorem 1]. Consider iteration (9) ( $\lambda=1$ ) with $n=m=2$ and design space $\mathcal{X}=\{x_{i}=(1,z_{i})^{\top},i=1,2\},z_{1}=-z_{2}=1$ . It is easy to show that, for any $w^{(t)}=(w_{1},w_{2})\in\Omega$ , iteration (9) maps $w^{(t)}$ to $w^{(t+1)}=(w_{2},w_{1})$ . Thus, unless $w_{1}=w_{2}=1/2$ to begin with, the algorithm alternates between two distinct points. This appears to be a rare example, as (9) usually converges in practical situations.

Further remarks and illustrations

One can think of several reasons for the wide interest in Algorithm I and its relatives. Similar to the EM algorithm, Algorithm I is simple, easy to implement and monotonically convergent for a large class of optimality criteria (although this was not proved in the present generality). Algorithm I is known to be slow sometimes. But it serves as a foundation upon which more effective variants can be built [see, e.g., Harman and Pronzato (HP07 2007) and Dette, Pepelyshev and Zhigljavsky (DPZ08 2008)]. While solving the conjectured monotonicity of (9) holds mathematical interest, our main contribution is a way of interpreting such algorithms as optimization on augmented spaces. This opens up new possibilities in constructing algorithms with the same desirable monotonic convergence properties.

As a numerical example, consider the logistic regression model

The expected Fisher information for $\theta$ from a unit assigned to $x_{i}$ is

We compute locally optimal designs with prior guess $\theta^{*}=(1,1)^{\top}(m=2)$ , and design spaces,

The design criteria considered are $\phi_{0}$ (for D-optimality) and $\phi_{-2}$ . We use Algorithm I with $\lambda=1$ , starting with equally weighted designs.

For $\phi_{0}$ , Corollary 1 guarantees monotonic convergence. This is illustrated by Figure 1, the first row, where $\phi_{0}=\log\det M(w)$ is plotted against iteration $t$ . Using the convergence criterion (4) with $\delta=0.0001$ , the number of iterations until convergence is 93 for $\mathcal{X}_{1}$ and 2121 for $\mathcal{X}_{2}$ . The actual locally D-optimal designs are $w_{1}=w_{20}=0.5$ for $\mathcal{X}_{1}$ and $w_{1}=w_{23}=0.5$ for $\mathcal{X}_{2}$ , as can be verified using the general equivalence theorem. This simple example serves to illustrate both the monotonicity of Algorithm I (when Theorem 1 applies) and its potential slow convergence.

For $\phi_{-2}$ , although Algorithm I can be implemented just as easily, Theorem 1 does not apply because the concavity condition (7) no longer holds. Indeed, Algorithm I (with $\lambda=1$ ) is not monotonic, as is evident from Figure 1, in the second row, where $\phi_{-2}=-\operatorname{tr}(M^{-2}(w))$ is plotted against iteration $t$ . This shows the potential danger of using Algorithm I when monotonicity is not guaranteed.

Although Theorem 1 does not cover the $\phi_{p}$ criterion for $p<-1$ , it is still possible that monotonicity holds for a smaller range of $\lambda$ . Calculations in special cases lead to the conjecture [Silvey, Titterington and Torsney (STT78 1978)] that Algorithm I is monotonic if $0<\lambda\leq 1/(1-p)$ . Theorem 1 provides further evidence for this conjecture, but new insights are needed to resolve it.

We have focused on local optimality. An alternative, Bayesian optimality [Chaloner and Larntz (CL89 1989), Chaloner and Verdinelli (CV95 1995)], seeks to maximize the expected value of $\phi(M(\theta;w))$ over a prior distribution $\pi(\theta)$ . The notation $M(\theta;w)$ emphasizes the dependence of the moment matrix on the parameter $\theta$ . It would be worthwhile to extend our strategy in Section 3 to Bayesian optimality, and we plan to report both theoretical and empirical evaluations of such extensions in future works.

Acknowledgment