Optimal Design of Experiments: Numerical Methods

W O R K I I V G P A P E R

OPTIMAL DESIGTJ OF EXPERIMENTS:

NUMERICAL METHODS

V.V. Fedorov

O c t o b e r 1986 WP-86-055

I n t e r n a t i o n a l I n s t i t u t e for Applied Systems Analysis

NOT FOR QUOTATION WITHOUT THE PERMISSION OF THE AUTHOR

Optimal Design of Experiments: Numerical Methods

KV; Fedorov

October 1 9 8 6 Cv'P-86-55

Working Papers are interim reports on work of the International Institute f o r Applied Systems Analysis and have received only Limited review. Views or opinions expressed herein do not necessarily r e p r e s e n t those of the Institute or of its National Member Organizations.

INTERNATIONAL INSTITUTE FOR APPLIED SYSTEMS ANALYSIS A-2361 Laxenburg, Austria

Acknowledgement

I am grateful to Prof. R . Munn for encouraging discussions and construc- tive remarks.

PREFACE

Optimal experimental designs became a r a t h e r efficient tool in applica- tions. T h e r e are numerous catalogues of optimal designs f o r some s t a n d a r d situations, f o r instance, when t h e r e s p o n s e function is a multidimensional power polynomial o r trigonometrical s e r i e s . In cases when t h e r e s p o n s e function is not of t h e approximation t y p e b u t i t s s t r u c t u r e is b a s e d on some

"physical" assumptions, o n e c a n not hope t o find a n optimal design suitable f o r specific conditions, a n d usually one needs t o a p p l y numerical methods t o find t h i s design.

This s h o r t s u r v e y is devoted t o numerical methods f o r t h e construction of optimal designs f o r e x p e r i m e n t s when a system "object u n d e r investigation

-

p r o c e s s of observation" is described b y t h e model (see (1) in t h e p a p e r ) which linearly d e p e n d s upon unknown p a r a m e t e r s a n d contains a n additional s t o c h a s t i c component (usually r e f e r r e d t o as a n e r r o r of observation, b u t i t could also r e f l e c t t h e s t o c h a s t i c n a t u r e of t h e object). The objective of a n e x p e r i m e n t is t o estimate t h e unknown p a r a m e t e r s a n d an optimal design (for instance, optimal location of observations) h a s t o provide t h e smallest e r r o r s of estimates. Usually. t h e s e e r r o r s are c h a r a c t e r i z e d b y t h e variance-covariance matrix (or some functions of it) which are t h e inverse of t h e so-called information matrix used in t h e paper.

Sections 1 a n d 2 contain some general information on numerical methods.

The subsequent s e c t i o n s deal with more specific situations. For instance, section 3 a n d 6 (containing some new r e s u l t s b y t h e author) deal with spatial- ly d i s t r i b u t e d observations a n d could be especially useful in optimization of monitoring systems, a v e r y a c u t e problem f o r many environmental studies. In Section 4, one c a n find algorithms which c a n be used f o r t h e design of e x p e r - iments r e l a t e d t o remote sensing of t h e e a r t h ' s atmosphere b y a satellite

r e

diometer (findings of optimal frequency bands o r "windows").

Prof. M. A. Antonovsky

CONTENTS

1

.

Introduction

...

¹

2

.

First Order Iterative Procedures

...

2 3

.

Construction of Optimal Designs Under Constraints

...

4 4

.

Optimal Designs when Controls Belong

t o a Functional Space

...

⁶

5

.

Discrete Designs

...

8 6

.

References

...

9

. vii

.

Opthal Design of Ekperiments: Numerical Methods

KV; Fedorov

In this paper numerical approaches for the construction of optimal designs will be considered for experiments described by the regression model

where f (x ) is a given set of basic functions, x EX , and X is compact; at least some of the variables x can be controlled by an experimenter, ^2p ERm are es- timated parameters, y* ER' is the i -th observation, and &* ER' is the ran- dom error, E [&*]

=

0, E ^[&* &,]

=

d t j . In practice, technically more compli- cated problems could be faced (for instance, y* could be a vector or errors could be correlated) but usually the methods are straightforward generaliza- tions of the methods developed for problem (1).

The most elegant theoretical results and algorithms were created for a continuous (or approximate) design problem when a design is considered to be a probabilistic measure defined on X , and an information matrix is defined by an integral M(C)

= jj

(x)*(x)€(dx). In this case, the optimal design of the experiment turns out to be the optimization problem in the space of pro- bability measures:

where @ is the objective- function defined by an experimenter.

The first ideas on numerical construction of optimal designs can be found in the pioneer works by Box and Hunter (1965) and Sokolov (1963), where some sequential designs were suggested. These procedures can be con- sidered as very particular cases of some iterative procedures for optimal design in construction, but nevertheless they implicitly contain the idea that one can get optimal design through improving intermediate designs by transferring a finite measure to some given point in X a t every step of the sequential design.

This idea was developed and clarified by many authors and the majority of algorithms presented in this survey (which does not pretend to be a his- torical one) are based on it.

2. First-rder iterative procedures.

I t will b e assumed t h a t

(a) t h e functions f ( z ) a r e continuous on compact X, (b) @ (M) is a convex function,

(c) t h e r e e x i s t s q such t h a t

and

(d) for any [EE(q) and any o t h e r

If t h e s e assumptions hold, t h e n t h e following iterative procedure will converge t o an optimal design:

where [(z,) is a design with t h e measure totally concentrated a t t h e point

Z s *

2,

=

Are m i n [ ~ (z:,

tS 1,

- v ( z ~ , t s ) l , (5)

Xs is t h e supporting set of t h e design [, , a, =7,, when z, =z,

+

^, and a,

=

^-%in [y, , pst / ( 1 -pSc)] , pst is a measure f o r point zsi of design

Ex.

To provide weak convergence, t h e sequence t y , j has t o satisfy, for in- stance, t h e following condition: 7, + 0 and C ys + m , s + m. In addition, some o t h e r alternatives for t h e sequence 17, j can b e found in Ermakov, 1983;

Fedorov, 1972; Fedorov and Uspensky, 1975; Fedorov, 1981; Silvey, 1982, Wu and Wynn, 1978. The iterative procedures (4), (5) comprise practically all t h e first-order methods widely discussed in t h e statistical literature since t h e late nineteen sixties. I t should b e pointed out t h a t t h e iterative pro- cedure can b e realized in practice if t h e optimization problem (5) is not very difficult from a computational point of view, i.e., if t h e dimension of X is not too high. I t is especially difficult to work with cases when t h e controllable variables belong to some functional space (see section IV).

There is a simple idea behind t h e iterative procedure (4), (5). If one wishes t o move along t h e "best" direction

t h e n for sufficiently small a, ( s e e (3)):

T h e r e f o r e

c a n b e chosen from t h e set of point measures t ( z ) , z EX a n d h a s t o be concen- t r a t e d a t t h e point:

The same idea is behind "deleting" some points from design

t,.

^These

points @,-j a r e "worst " in t h e s e n s e (6).

For fulfillment of t h i s f a c t , assumption (d) is crucial. The majority of op- timality criteria used in p r a c t i c e s a t i s f y t h i s assumption. But f o r some quite natural c r i t e r i a , f o r i n s t a n c e ,

where z , is given.

Q

is t h e variance of

G T f

(zJ a n d '"" s t a n d s f o r pseu- doinversion, formula (3) is not generally valid. One c a n s t i l l use t h e i t e r a t i v e p r o c e d u r e (4), (5) applying t o some regularized version of t h e initial problem:

where M ( t J is regular matrix (M ([J

a).

Then

To adjust t h e i t e r a t i v e p r o c e d u r e (4), (5) t o p a r t i c u l a r optimality c r i t e r i a , t h e following formulae c a n be useful:

where t h e e x i s t e n c e of a corresponding derivative is assumed,

The convergency rate of the above mentioned algorithms decreases in the vicinity of the optimum. As in general optimization theory, attempts were made to develop second order methods. These methods are based on quadrat- ic approximations of the function @ [ M I and it is necessary to assume the ex- istence of derivatives

am/

apt

.aZ@/

apt

ap, .

where p t is the weight for sup- porting point x t . Second-order algorithms have at least two features which are handicaps for their use in practice: first, at every step it is necessary to invert the matrix

I$@/

^3pt

apt 1,

and second, all existing modifications can handle only the discrete operability region X (see Ermakov (ed), 1983 Ch 4, Wu, 1978).

In the late seventies, some attention was paid to algorithms which work in the space of information matrices M (X); they are computationally effective if one can easily find the mapping E(X)

+M

(X) (see, for example, Gribik and Kortanek, 1977). B u t usually it is very difficult to realize this mapping nu- merically.

3. Construction of _optim a1 designs under constraints.

In (2) there is only one constraint j [ ( d x )

=

1. If one considers addi- tional constraints, say j + ( x ) [ (dx)

* .

^{where +(z)} is a vector of given functions (c

,+ ^mk),

then the iterative procedure becomes technically more complicated, although based on the same ideas, see (Fedorov and Gaivorov- sky, 1984). In this case, (6) became equivalent to the following optimization problem

subject to j + ( z ) i ( d x ) * . X

Due to the classical theorem of Caratheodory, design

Ts

can be found within the set of designs containing no more than k +1 supporting points.

The optimization problem (8) is essentially more complicated than (5). From a computational point of view, the dual problem (see Karlin and Studden, 1966, ch. XII):

where porting points of U = f u :ut 10-,i = l , k

t, .

They have to coincide with the solutions x i , x z j, can be useful to define the location of the sup-

. ,...

^of

(9). The corresponding measures can be found by any linear programming al- gorithm.

In some applications (see Wynn, 1982), designs have to be restricted in the following sense

f

x

1 J 1

for all A ~ 1 ( (10)

A A

where

f

to(dz)=c, c 21, A C1( and a measure [, is atomless (see Karlin and X

Studden, 1966, p. 233). In this case (6) leads t o t h e following problem

where a probability measure

js

has t o satisfy (10).

I t is evident t h a t

7,

has t o coincide with

to

on any s u b s e t of X, where

~ ( x , [ , ) a and has t o b e equal 0 otherwise (compare with theorem 1 from Wynn, 1982). Computationally, t h e s e a r c h f o r t h e s e sets c a n b e realized through discretization of X . The idea of t h e iterative procedure (4). (5) will apply once again if one will additionally delete from design

e,,

^{some sets}

where ~ ( x ,c,)>O. Thus

where Es is t h e set of new included points and Ds is t h e set of deleted points.

The procedure similar t o (12) (but without t h e operation of deleting) was considered by Gaivorovsky (1985), and i t s weak convergency was proven under r a t h e r mild conditions. Usually, deleting "bad" points essentially im- proves t h e quality of t h e iterative procedures (compare with t h e traditional case, Atwood, 1973, Fedorov and Uspensky, 1975).

Let E(cJ b e a class of probability measures

c

with supporting s e t s A

cX,

and [(dz)=c,(&), when x EA and equal t o 0 otherwise. For any design prob- lem (2), (lo), t h e r e exists a n optimal design c'E$(cJ, see, f o r instance.

Wynn, 1982, whose results have t h e i r origin in t h e classical moment space theory, particularly in t h e Liapunov Theorem, see Karlin and Studden, 1966, C h . . VIII). Therefore, i t is reasonable t o demand t h a t

cs

=c(A, ,dx ) €Z(c,J,s =1,2,

...

Iterative procedure (12) does not satisfy t h e l a t t e r demand. Instead of t h i s t y p e of iterative procedure (which r e p e a t s t h e idea of (4), (5)), one c a n apply a n iterative procedure of t h e exchange type:

where

~ C , ( ~ ) = ~ F , ( ~ X ) = ~ ~ . As

M s

=O.

E8 0 8

Ds

a,,

E,

nD,

=O.

If

(e) derivatives

&

e x i s t ,

(f)

to

has a continuous density p,(x),

and assumptions (a)-(c) hold, t h e n f o r sufficiently small 6,. t h e approxi- ma tion

@[M (cs +I)] =@[M(

cs

)I +[7(x,+. C,) -7(xs

- . ^ts

^{)I 6, +o (6,)- (15)}

where X$EE,, x;€LlS and y(x

.o=f

~ ( x ) i ( c ) f ( x ) can b e used. If X is

covered b y t h e g r i d Xd with d e n s i t y proportional t o p,(z), t h e n

12

^{a n d xs}

-

c a n coincide with i t s nodes a n d Es , Ds with some cells of t h i s grid.

From (15), it is c l e a r t h a t t o provide approximately t h e s t e e p e s t d e s c e n t on e v e r y s t e p of t h e discretized version of p r o c e d u r e (13). one h a s t o find

x:

AT^

min y ( z ,

ts

⁾ a n d zs-=Arg max y ( z ,

ts

),

xEX6\A6r ad*

where A ds is a discrete analogue of A,. I t is worthwhile t o point out t h a t f o r t h e discretized version of t h e i t e r a t i v e p r o c e d u r e , one c a n use a recursive formula f o r M , - . ~ (see (17)) t o simplify calculations. Complementing assump- tions (a)-(f) b y assumptions:

(g) f o r a n y design

t

with Q[M (t)]SQ ^<m a n d a n y C t,IA:y(x,t)=C1=0

.

(h) 6s 4,

C

6, ^+m

S

t h e weak convergency:

lim @[M(ts)]= min @[M(t)]

s ^-+m €cE(t0)

c a n b e proven.

Result (17) is based on t h e f a c t t h a t t h e fulfillment of t h e inequality:

max y(x ,[*)( min y ( z , t * ) z

a*

x

a\x*

is a necessary a n d sufficient condition f o r a design

t*

t o b e optimal (x* is a supporting set of

t*

^).

4. Optimal designs when controls belong to a functional space.

This case will b e illuminated h e r e b y a r a t h e r specific example which nevertheless r e f l e c t s t h e major difficulties.

Let f (v)ERm a n d x = j f ( v ) h ( v ) d v , where h ( v ) c a n b e controlled.

v

h (v) EH. If one manages t o c o n s t r u c t t h e mapping X ( H ) cRrn of t h e set H, t h e n all a p p r o a c h e s discussed in t h e previous sections c a n b e used without a n y a l t e r a t i o n s t o find a n optimal design on X. The problem t o b e f a c e d a f t e r w a r d s is t h e construction of a n inverse mapping X+H t o c o n v e r t

t i

t o some design defined on H. The l a t t e r problem is beyond t h e s c o p e of t h i s p a p e r a n d its discussion c a n be find in Kozlov, 1981, Ermakov (ed), 1983, Ch.7. For t h e c a s e discussed in section 2 , t h e situation is slightly simpler be- c a u s e of t h e Equivalence Theorem (see, f o r instance, Ermakov (ed), 1983 Ch.2); t h u s only boundary points X(H) of X(H) are needed f o r optimal design construction. Unfortunately, t h e numerical construction of X(H) h a p p e n s t o b e sometimes a very difficult problem a n d it could b e more efficient t o work in t h e original s p a c e H.

To be more specific, let us assume that VEVCR' and OSh ( v ) S and res- trict ourself to the design problem in section 2. The most straightforward approach consists of discretization of V and approximation of h (v) by same piecewise function. Under rather mild conditions, it is possible to prove that there exists optimal design which supporting points belong to the set H=lh (v):h (v)[l -h (v)]=O,vEVj, see, for instance Fedorov, 1806. If V is discretized by a grid with elements A,, then the simplest version of pro- cedure (4), (5) (without "deleting" operation) converts to the following one:

a)

ts

+1=(1 -as

>t,

+a, t(hs) ;

b) steps for finding h, :

-collect all A j which negatively contribute to the sum

where Fj

=I

^f (v)dv (usually Fj sf (vj)Aj).

* j

-

put h, ( v ) = l , vEAj if A j was chosen in the previous stage, otherwise hs(v) =O ,

-

the fulfillment of the inequality p(C,)<t+

q(ts)M

(c,) tests that h, can be used for C(h, ).

This iterative procedure guarantees that

where

t i

is an optimal design for the discretized design problem and can be called a A

-

optimal design. When a c v s b and functions f (v) constitute a Tchebycheff system over the open interval (a,b), where a and b are possibly infinite, then the rather effective iterative procedure can be used for op- timal design construction. The idea of this procedure is based on the follow- ing result (see, for instance, Fedorov, 1986).

Let h

(v)EH

and let I be the number of separate nondegenerate intervals where h ( v ) = l with the special convention that an interval whose closure contains point a or b, is counted as 1/2. For any point x =,I* stands for the least possible I. Then a necessary and sufficient condition that z belongs to the boundary of X is that I* ((m -I)/ 2. Moreover, every boundary point corresponds to a unique h (v ) with I (z )=I * ( x ).

Let now

G

=(v l,...,vm where a S v l S

. .

Sum -lSb. According to the previous result, there exist optimal designs with all supporting points (in the operability region H ) which have the following structures:

h ( v ) = 1 . ~ ~ ( a , ~ ~ ) : 0 . ~ ~ ( ~ , , v ~ ) ; 1 , v ~ ( v ~ , v 3 > ;

. . -

and _h ( v ) =1 -h (v)

.

That fact allows for modification of the iterative procedure (4), (5) (without deleting "bad" points) to the procedure with maximization in space.

with dimension less than or equal to (m -I), where m is a number of basic

functions:

h, =~rgmi_n@[x 7(%),ts

I,

Y*V

b b

-=b .xi(%)

= Sf

(v) h ( u ) d v and x Z ( ~ ) = S f ( v ) b (v)dv where a <vlS .

. .

^<V,

a a

The design problem considered in this section comprises t h e major diffi- culties which can b e m e t when X is a functional space. Other examples can b e found in Mehra, 1976, Pazman, 1986. These authors use mainly t h e same ideas surveyed h e r e . In concluding this section, i t would b e worthwhile t o notice t h a t parametrization of controls (a r a t h e r standard method in optimal control theory), e.g. linear approximation cSTq (v) of h (v) in our example, could b e a useful tool allowing one t o convert t h e original problem t o a finite dimension design problem.

5. Discrete designs

To c o n s t r u c t optimal d i s c r e t e (or e x a c t ) designs, a number of exchange t y p e algorithms can b e used (for detailed information, see Cook and Nachtsheim, 1980; Johnson and Nachtsheim, 1983; Steinberg and Hunter, 1984).

The idea of t h e simplest algorithm (originated b y Mitchell, 1974) can b e formulated in t h e following way:

A f t e r t h e s-th s t e p t h e r e is a design

tNs

=)xis,

...,

xNs

1

, where some sup- porting points c a n coincide. This design is complemented by k points:

Then t h e same number of points:

a r e deleted and one arrives a t t h e new design

IN,,

containing N observa- tions. The notation t K + x ( o r z) means t h a t a point is included in (or ex- cluded from) design tK,Xs comprises all of t h e different supporting points of t h e design from t h e previous stage.

In practice, t h e excursion length k is usually r a t h e r modest (1-3) and t h e r e a r e no indications t h a t a n increase could b e useful. Iterative pro- c e d u r e ( l a ) , (19) are computationally simple and often lead t o very good results, especially when one f a c e s d i s c r e t e X , for example,

z

,=*I.

In t h e iterative procedures (18), (19). t h e deletion and complementary s t e p s a r e separated. If w e unite them (Fedorov, 1972), t h e n w e arrive a t t h e following iterative procedure (with excursion length 1):

where

Xs

is t h e supporting set of

CN,.

This procedure demands

N /

2 times more calculations a t every s t e p than (18), (19), but in most cases i t gives b e t t e r final results, see Johnson and Nachtsheim, 1983. The above minimiza- tion problem is equivalent t o coordinate wise minimization of

@[MI:

minmin 9[M(zi,

,...,

xjs

,...

z ~ , ] f ^{z j a}

if one s t a r t s t h e numerical optimization in (21) with zP=zj,.

A similar choice of a n initial point is a p p r o p r i a t e in many optimization problems, b u t not in t h e optimal design of experiments when t h e objective function usually has a large number of local minima along t h e variation of z j ^, and (21) will lead t o t h e local minimum closest t o xjs. The application of (20) helps t o approach t h e global minimum b y explicit forcing of z ⁺ t o b e away of x j-. Procedures (18), (19) o r (20) become a practical tool when one manage t o find a simple formula f o r calculation of increments f o r

@[MI

at every stage.

For t h e majority of widely used c r i t e r i a (D-criterion, linear c r i t e r i a and s o on) t h e s e formulas can b e found in t h e above c i t e d publications (see also sec- tion 2).

In s p i t e of t h e r a t h e r long history of t h e numerical procedures dis- cussed above, t h e i r convergence properties a r e not well known e x c e p t f o r numerous empirical results. I t is not a problem, f o r instance, t o prove t h e convergence of (20) t o some design b e t t e r than a n initial one b u t i t has not y e t been proven t h a t t h e limit design has t o b e optimal.

Atwood, C.L. (1973) Sequences Converging t o D-Optimal Designs of Experi- ments, Ann.Stat., Vol.1:342-352.

Box, G.E.P. and Hunter, W (1985) Sequential Design of Experiments f o r Non- linear Models. Proc. IBM Scientific Computing Symp. Statis tics, 113-120.

Cook, R.D. and Nachtsheim, C.J. (1980) A Comparison of Algorithms for Con- s t r u c t i n g Exact D-Optimal Designs, Technometrics, Vo1.22:682488.

Ermakov, S.M. (ed.) (1983) Mathematical Theory of t h e Design of Experiments (in Russian), p.386.

Fedorov, V.V. (1972) Theory of Optimal Experiments, Academic Press, N.Y., p.

292.

Fedorov, V.V. and Uspensky, A.B. (1975) Numerical Aspects of Design and Analysis of Experiments, Moscow S t a t e University, Moscow (in Russian), p. 167.

Fedorov, V.V. (1981) Active Regression Experiments in Mathematical Methods of Experimental Design, ed. Penenko, V.V., Nauka (in Russian), p.19-73.

Gaivoronsky, A. and Fedorov V. (1984) Design of Experiments under Con- s t r a i n t s , WP-84-0, IIASA, Laxenburg, Austria, p. 11.

Gaivoronsky, A. (1985) Stochastic Optimization Techniques f o r Finding Op- timal Sub-measures, WP-85-28, IIASA, Laxenburg, Austria, p.52.

Gribik, P.R. and Kortanek, K.O. (1977) Equivalence Theorems a n d Cutting Plane Algorithms f o r a Class of Experimental Design Problems, SIAM J.

Appl. Math., Vo132:232-259.

Johnson, M.E. a n d Nachtsheim, C.J. (1983) Some Guidelines f o r Constructing Exact D-Optimal Designs on Convex Design Spaces, Technometrics, Vo1.25:271-277.

Karlin, S. and Studden, W.J. (1966) Tchebycheff Systems: with applications in Analysis a n d Statistics, Wiley a n d Sons, N.Y., p. 578.

Kozlov, V.P. (1981) Design of Regression Experiments in Functional Spaces, in Mathematical Methods of Experimental Design, ed. Penenko, V.V., Nauka (in Russian), p.74-101.

Mehra, R.K. (1976) Synthesis of Optimal Inputs f o r Multiinput- Multioutput Systems with Process Noise, in "System Identification: Advances and Case Studies" ed. D.K. Mehra, D.G. Lainiotis, Academic P r e s s , N.Y. p p . 211-250.

Mitchell, T. J. (1974) An Algorithm f o r t h e Construction of D-Optimal Experi- mental Designs, Technometrics, Vo116:203-210

Pazman. A. (1986) Foundations of Optimum Experimental Design, VEDA, D.

Reidel Publishing Company, p.228.

Silvey, S.D. (1980) Optimal Design, Chapman a n d Hall, London, p.86.

Steinberg, D.M. a n d Hunter W.G. (1984) Experimental Design: Review a n d Com- ment, Technometrics, Vo1.26:71-97.

Sokolov, S.N. (1963) Continuous Design of Regression Experiments. Teor.

Verejatuost.; Primen., Vo1.0:95-101 (a), 318-323 (b).

Wu, C.F. (1978) Some Iterative Procedures For Generating Nonsingular Optimal Designs, Commun. S t a t i s t . (Theor. Math.), Vol. A7(14):1399-1412.

Wu, C.F. and Wynn, H. (1978) The Convergence of General Step-Length Algo- rithms f o r Regular Optimum Design Criteria, Ann. Statist., Vol. 6:1273- 1285.

Wynn, H. (1982) Optimum Submeasures with Applications t o Finite Population Sampling, in "Statistical Decision Theory and Related Topics 111", Vo1.2, Academic P r e s s , N.Y., pp. 485-495.