International Journal of Scientific & Engineering Research, Volume 6, Issue 1, January-2015 1847

ISSN 2229-5518

Reference Number: 41127004; Paper ID: I058895

SEGMENTED LINE SEARCH METHOD FOR CONSTRUCTING

D-OPTIMAL EXACT DESIGNS IN CONTINUOUS EXPERIMENTAL AREAS

F.I Nsude

Department of Statistics

Akanu Ibiam Federal Polytechnic, Unwana, Ebonyi State, Nigeria

P.E. Chigbu

Department of Statistics, University of Nigeria, Nsukka, Enugu State, Nigeria

Abstract

An effective segmented line search technique for constructing D-optimal exact designs in a polynomial response
function of degree z is presented in this paper. The direction of search is a weighted average from all segments of the experimental trials; the weight being proportional to the mean square errors from the segments. The method is rapidly convergent and its application in continuous experimental areas is shown.

Key words: variance function, Exact Designs, Continuous design, Exchange algorithm.

~

1. Introduction

The method reported here is a segmented line search
procedure for constructing D-optimal exact designs in a continuous experimental region. In this search system, variance of polynomial response function of degree z is approximated by a polynomial function of degree 2z. Maximizer of this polynomial is determined by segmenting the experimental area into segments. Thereafter, this maximizer is used to replace the point of minimum prediction variance amongst the points

X is a collection of factor levels x , which constitutes a compact, continuous and metric space, and at each

point, 𝑥 ∈ 𝑋� , the response is observed with a random
error, 𝑒𝑥 , which is a point in the space, ∑ x which is a
non-negative continuous function, define in 𝑋� (Chigbu
and Oladugba, 2010) [4].
In this study, treatment or support point, x , was
already in the design. The search method is via the
selected such that each element of

xi , is attached a

Super Convergent line Search Series reported by
Onukogu (1997) [8] which has shown to have rapid convergence capability; see Onukogu and Chigbu (2002, p. 112) [9], and Ugbe and Chigbu (2014) [11].
Every line search technique has four characteristics
weigh wi; wi ≥ 0 and ∑n 1 wi = 1, where n is the

i=

number of support points in the design. From a
collection of n support points we form a design matrix,

X, as well as normalized Fisher’s information matrix,

𝑋 ′ 𝑋 𝜎2

sequence as follows: (a) the starting point, 𝑥̅, which is
M( ξ ) =

𝜉 𝜉 , where X is the design matrix of the

𝑛

an n-component vector, x ; (b) an n-component

model terms (the columns) evaluated at specific

′ 𝑋𝜉 is

direction vector, 𝑑; (c) the step-length,

ρ; and (d)

point/treatment in the design space (the rows), 𝑋𝜉
Movement to the iterate, ( 𝑥 = 𝑥̅ + 𝑑𝜌 ). This
a non-singular matrix (information matrix) and

′ −1

set of sequential activities applies as well to this
𝜎 2 (𝑋𝜉 𝑋𝜉 )
is the variance-covariance matrix of the
technique herein reported.
least square estimate of the response parameters.
This search method also makes use of variance exchange of points procedure which is originally
Our interest is to find n-points design, 𝑥𝑖
∈ 𝑋� , 𝑖 =
reported by Fedorov (1972) [5], and modified severally over decades before the Variance Exchange Method reported by Atkinson and Donev (1989) [1].

2. Design Background Problem

1, 2, … , 𝑛, such that the resultant information matrix
(𝑋′ 𝑋𝜉 ) is maximized. A design criterion that
maximizes the determinant of the information matrix
|𝑋′ 𝑋𝜉 | or, equivalently, minimizes |𝑋′ 𝑋𝜉 |−1 is referred

𝜉 𝜉

~
Given the basic experimental triple X,
Fx , ∑ x },
to as D-optimality. Other optimality criteria have been
reported severally (Atkinson and Donev, (1992) [2]; Onukogu and Chigbu (2002) [9]).
where 𝐹𝑥 is the response functions, which is a space of
finite dimensional linear function and continuous in ~

IJSER © 2015 http://www.ijser.org

International Journal of Scientific & Engineering Research, Volume 6, Issue 1, January-2015 1848

ISSN 2229-5518

Therefore in this paper, the intention is to present a search method which seeks iteratively an n-point design
measure, 𝜉𝑛 that maximizes the determinant of M(𝜉∗ );
by 𝑐𝑏 . Matrix of mean square error of each of the kth
segment is given by
i.e Max| M(𝜉∗ )|; M(𝜉𝑛 ) ∈ 𝑆𝑝 × 𝑝 , where 𝑆𝑝 × 𝑝 is a set

−1 −1 −1

𝑀 = 𝑋 + 𝑋 𝑋′ 𝑋𝑏 𝑐 𝑐′ 𝑋′ 𝑏 𝑋 𝑋 (3.4)
of all non-singular p x p information matrices defined

𝑘 𝑘

𝑘 𝑘

𝑘 𝑏 𝑏

𝑘 𝑘 𝑘

in 𝑋� .
where 𝑋

= 𝑋′ 𝑋 , and 𝑀
is an (m+1) x (m + 1)

𝑘 𝜅 𝜅 𝑘

3. Theoretical Framework
The method is summarized by the following sequence of steps.
matrix, m being the number of factors/variables.

3.8. The average mean square error is minimized to

obtain the matrix of convex combination, 𝐻𝑘 ; Thus,

min (𝑀(𝑑𝑡)) = min(∑𝑠
ℎ𝑡𝑘 𝛼𝑡𝑘
) and as the cross

3.1. Pick at random an initial n-point design of

nonsingular information matrix.

3.2. Obtain its extended design matrix, X, according to the model terms, and the Fisher’s information matrix,

terms are zero; we have


min �𝑀(𝑑𝑡)� = 𝑚𝑖𝑛 ∑𝑠 ℎ2 𝑀(𝑎 ) (3.5)

𝑘=1 𝑡𝑘 𝑡𝑘

𝑋′ 𝑋 = 𝑀(𝜉).
where 𝑑𝑡 = ∑𝑠
0, 1, … , 𝑚, ℎ𝑡𝑘 ≥ 0.
ℎ𝑡𝑘 𝑎𝑡𝑘

𝑠

𝑘=1

ℎ𝑡𝑘
= 1; 𝑡 =

3.3. From the above, compute standardized generalized variance function of degree 2z as

𝑑�𝑥, 𝜉� = 𝑛𝑥′ 𝑀−1 (𝜉)𝑥, (3.1)
= 𝑔(𝑥) (3.2)

The values, ℎ𝑡𝑘 , are obtained by taking partial derivatives of 𝑀(𝑑𝑡 ) with respect to ℎ𝑡𝑟 of equation
(3.5), equating to zero, and solving.

3.9. The average information matrix from all the segments is given by

= 𝛽0 + ∑𝑚

𝛽𝑖 𝑥 + ∑𝑚

𝛽 𝑥

+ ∑𝑚−1 ∑𝑚

𝛽 𝑥 𝑥 …

′

1×𝑝

𝑖=1 𝑖

𝑖=1 𝑖 𝑖

𝑖=1

𝑖<𝑗 𝑖

2

𝑖 𝑗

2

𝑀𝑑 = 𝐻𝑋′ 𝑋 𝐻′

𝑠

′ ′

= (1 𝑥1𝑖 𝑥2𝑖 . . . 𝑥1𝑖 𝑥2𝑖 𝑥1𝑖 𝑥3𝑖 . . . 𝑥1𝑖
𝑥2𝑖 . . . )
= � 𝐻𝑘 𝑋𝑘 𝑋𝑘 𝐻𝑘

𝑘

Note that the number of factor interaction is 𝑚(𝑚−1),

2

where m is the number of factors or variables.
Where H is a matrix of convex combination, and H =
(H , H , . . . , H )

1 2 s

3.4. Partition the experimental region,𝑋� , either into

𝑋𝜄 𝑋 = 𝑑𝑖𝑎𝑔 (𝑋′ 𝑋 , 𝑋′ 𝑋 , … , 𝑋′ 𝑋 );

𝑆𝑘 ; 𝑘 = 1,2. . . , 𝑠 non-overlapping segments such that

1 1 2 2

𝑠 𝑠

𝑆1 ∩ 𝑆2 . . .∩ 𝑆𝑠 = ∅, and 𝑆1 ∪ 𝑆2 . . .∪ 𝑆𝑠 ≤ 𝑋� ; or with common boundary such that 𝑆1 ∩ 𝑆2 . . .∩ 𝑆𝑠 ≠
∅, and 𝑆1 ∪ 𝑆2 . . .∪ 𝑆8 ≤ 𝑋� : See (Ugbe and Chigbu
(2014) [11].
3.5. From each kth segment, pick or select 𝑛𝑘 supports

3.10. The response vector, z, is obtained from the variance function, g (x). Hence,
𝑧 = (𝑧0 , 𝑧1 , … , 𝑧𝑚 )′ ; 𝑧𝑖 = 𝑔�𝜂𝑖𝑗 �, (3.6)
𝑖 = 0, 1, … , 𝑚; 𝑗 = 1, 2, … , 𝑚;
points, 𝑛𝑘 ≥ 𝑚 + 1, 𝑘 = 1,2, . . . 𝑠; ∑𝑠
𝑛𝑘 = 𝑁0, the
total number of selected points in all the segments.
3.6. Obtain the design matrix, 𝑋𝑘 , 𝑘 = 1, 2, . . . , 𝑠, for each segment, and the information matrices,

′

Where 𝜂𝑖𝑗 is the ith and jth element of the average information matrix, 𝑀𝑑.

3.11. Obtain the vector, starting point (𝑥∗ ),

𝑗

𝑋 = 𝑋𝑘 𝑋𝑘 (3.3)

𝑁0 −1

∗ 𝑖

3.7. Compute 𝑋𝑏𝑘, the matrices of

𝑥 = � 𝑤𝑖 𝑥𝑖 ; 𝑤𝑖 = ∑𝑁0

;
𝑎−1
interaction/biasing effect of the variable for each
segment; where k = 1, 2, . . ., s; and the vector of
𝑎𝑖 = 𝑥𝑖

𝑖=1

𝑁 𝑖

𝑖=1 𝑖

coefficient of interaction/biasing terms of g(𝑥) is given

′ 𝑀−1 (𝜉

)𝑥 ; (3.7)
And the direction vector (𝑑∗ ) is,

IJSER © 2015 http://www.ijser.org

International Journal of Scientific & Engineering Research, Volume 6, Issue 1, January-2015 1849

ISSN 2229-5518

𝑑̂ = 𝑀−1

𝑧 (3.8)

∞

 �𝑥 �

𝑗 =1

= 𝑥𝑟
𝑑 = (𝑑0 , 𝑑1 , . . . , 𝑑𝑚 )′ and normalized such that

1

𝑑̂ ′∗ 𝑑̂ ∗ = 1. Thus, 𝑑∗ = 𝑑(𝑑′ 𝑑)−2; and the optimal step-
Proof


length (𝜌∗ ), 𝑔(𝜌) = 𝑔(𝑥

𝑗

+ 𝜌𝑑), and by taking
Define the information matrix at the jth iteration as
derivatives of the function with respect to 𝜌, we have

𝑀𝑗 (𝜉𝑛 ) = ∑𝑖=1 𝑥𝑖 𝑥𝑖 (𝜉𝑛 ): See Gaffke and Krafft(1982)

𝑑[𝑔(𝜌)]
= 0;

𝑛 ′

[7] , then
𝑑𝜌

𝑀𝑗 = 𝑥1 𝑥1 + 𝑥2 𝑥2 +. . . +𝑥𝑟 𝑥𝑟 +. . . 𝑥𝑛 𝑥𝑛

′ ′ ′ ′

3.12. The new point is 𝑥𝑗 = 𝑥

+ 𝜌∗ 𝑑∗ , 𝑗 = 1, 2, . . .

 ∑𝑛−1 𝑥 ′ + 𝑥 𝑥′

𝑖=1

𝑖 𝑥1

𝑟 𝑟

3.13. Take the point, 𝑥1, to step (3.5) above and add it
to the segment it falls into, so that the number of
support points for the segment is (𝑛𝑘 + 1), and follow the steps down to obtain a second point 𝑥2. Again, add the points, 𝑥2, to the segments it falls, and continue with the steps to obtain a third pint, 𝑥3 , etc. This process is continued until it converges to a point, 𝑥 ∗.
Set 𝑀0 = ∑𝑛−1 𝑥𝑖 𝑥 ′

𝑀𝑗 = 𝑀0 + 𝑥𝑟 𝑥′ ,

Then the determinant becomes

|𝑀𝑗 | = |𝑀0 + 𝑥𝑟 𝑥′ |

′ 𝑀−1 𝑥 )| (4.1)

3.14. Compute the prediction variance of each design

 |𝑀0 (1 + 𝑥𝑟 0 𝑟
points as 𝑥 ′ 𝑀−1
(𝜉0 )𝑥𝑖 ; i = 1, 2, . . . , n and select the
The information matrix at the (j-1) iteration is 𝑀 =
one with minimum prediction variance 𝑉�𝑥𝑚𝑖𝑛 � = ′

𝑗−1

′

𝑚𝑖𝑛

𝑀−1 (𝜉0 )𝑥𝑚𝑖𝑛
. Add the point 𝑥 ′ row vector to

𝑀0 + 𝑥𝑘 𝑥𝑘

augment the n rows extended matrix, and remove the
point with minimum predicted variance in the design to

 |𝑀𝑗−1 | = |𝑀0 + 𝑥𝑘 𝑥𝑘 |

have the determinant as
�𝑀𝑗−1 � = |𝑀0�1 + 𝑥𝑘

1 𝑥 �| (4.2)

′ 𝑀−

′ ′

�𝑀(𝜉0 ) + 𝑥𝑗 𝑥𝑗 − 𝑥𝑚𝑖𝑛 𝑥𝑚𝑖𝑛 �

= |𝑀(𝜉0 )|�1 + ∆�𝑥𝑗 , 𝑥𝑚𝑖𝑛 ��
By comparing equation (4.1) and (4.2), the sequence moves from a point of low variance to a point of high variance, therefore
Note that 𝑥𝑚𝑖𝑛 is the point of minimum prediction

𝑥𝑘

0 𝑘 𝑟

0 𝑟

variance in the design and 𝑥 ∗
is the point of maximum

′ 𝑀−1 𝑥

≤ 𝑥′ 𝑀−1 𝑥

′ 𝑀−1 𝑥

≤ 𝑥 ′ 𝑀−1 𝑥 , (4.3)
prediction variance from the experimental region, 𝑋� .

 𝑥𝑗−1

𝑗−1

𝑗−1

𝑗 𝑗 𝑗

3.15. With the new point in the design, compute the

Since 𝑀0 = ∑𝑛−1 𝑥𝑖 𝑥′ is the same in both equations
inverse information matrix, and the variance function,
�𝑀0 �1 + 𝑥𝑟

𝑖=1 𝑖

0 𝑟 0

𝑘

0 𝑘

𝑔�𝑥�, and continue the other processes of the steps to
obtain another converging point, which is exchanged
with a point of minimum prediction variance from the

′ 𝑀−1 𝑥 �� ≥ �𝑀 �1 + 𝑥′ 𝑀−1 𝑥 ��

 |𝑀𝑗 | ≥ |𝑀𝑗−1 |

′ −1

′ −1

design.
Then, since 𝑥𝑘 𝑀0

𝑥𝑘 ≤ 𝑥𝑟 𝑀0

𝑥𝑟 , it follows from

′ −1 ∞

equation (4.3) that the point {𝑥𝑗 𝑀𝑗

𝑥𝑗 }𝑗=1, is

3.16. Is �det 𝑀𝑗 (𝜉)� − det 𝑀𝑗−1 (𝜉) |
≤ 𝜖; 𝜖 ≥ 0?
Yes; stop and 𝑀𝑗 (𝜉) is D-optimal.
4. Convergence of the sequence
monotonically increasing, and is sure to converge to a
point, 𝑥𝑟 .
=> {𝑥𝑗 }∞ = 𝑥 . Therefore, the Kolmogrov’s condition
for convergence with probability one is assured, see
Feller (1966, pg 259), and Onukogu and Nsude (2014) [10].

Given a line sequence 𝑥 = 𝑥

𝑗

+ 𝜌𝑗 𝑑𝑗 ; 𝑗 = 1, 2, … ;

The sequence has the following attributes:

∞ ∞

then, �𝑥 �

𝑗=1

= {𝑥

𝑗

+ 𝜌𝑗 𝑑𝑗 }𝑗=1

IJSER © 2015 http://www.ijser.org

International Journal of Scientific & Engineering Research, Volume 6, Issue 1, January-2015 1850

ISSN 2229-5518

(i) The direction vector is a weighted average from all segments; the weight being proportional to the mean square errors from the segments; see equation (3.7) and the sequence converges with minimax point; see Atkinson and Donev (1992, chapt. 11) [2].
(ii) Again, the response function was not use, only its variance function.
(iii) Division of experimental region into s over-lapping or non-overlapping segments avoids the use of second derivative in obtaining Hessian matrix which at times constitutes problem of inevitability in
Thus, the variance function is
𝑔�𝑥� = 2.5624 − 3.8249𝑥2 + 5.0624𝑥4
Now partition the experimental region −1 ≤ 𝑥 ≤ 1 into
k=2 non over-lapping segments; −1 ≤ 𝑆1 ≤ −0.25,
and 0.25 ≤ 𝑆2 ≤ 1 and in the kth segment, 𝑔(𝑥) is
approximated by a first-order linear function.
Thus, 𝑦𝑘 = 𝛽00 + 𝛽0𝑘 𝑥 ; 𝑘 = 1, 2.
Pick points in the 𝑠𝑘 segments, and obtain the 𝑋𝑘 design
matrices as thus:
most line search algorithms, eg. Newton’s method is achieved.

5. Numerical Demonstrations

1

𝑋1= �1

1

1

−1.0000

−0.7500

−0.5000

−0.2500

1

� ; 𝑋 2 = �1

1

0.2500

0.5000 �

0.7500

1.0000

The problems below are used to illustrate the numerical
and their information matrices are respectively given as
behavior of the search system.
𝑋�
4 −2.5000

Problem 1: Consider a polynomial regression function,

1 = �−2.5000 1.8750 � ;

2 = � 4 2.5000

𝑓�𝑥� = 𝛽0 + 𝛽1 𝑥 + 𝛽2 𝑥2 , −1 < 𝑥 < 1; obtain a 4-
point exact D-optimal design. (Atkinson and Donev,
𝑋�
2.5000 1. � ;
1992) [2]; (Atkinson, Donev and Tobias, 2006) [3] used
where 𝑋 �

′ 𝑋

and 𝑋�

′ 𝑋

this problem to find exact D-optimal design of a

1 = 𝑋1 1

2 = 𝑋2 2

quadratic polynomial of a single factor. Begin with
initial design points as
𝜉0 = {−1, −0.3333, 0.3333, 1}
The matrices of biasing terms, 𝑥2 , 𝑥4 , of the variance
function, 𝑔�𝑥�, of each segment are
The extended design matrix is,
1 −1 1
𝑋(𝜉0 ) = �1 −0.3333 . 1111� ;
𝑋𝑏1
⎛1.0000
= ⎜0.5625
0.2500
⎝0.0625
1.0000⎞
0.3164⎟ ;
0.0625
0.0156⎠
1 0.3333
1 1
. 1111
1
⎛0.0625
0.0156⎞
𝑋′ 𝑋(𝜉0 ) = 𝑀(𝜉0 )
4.0000 0 2.2222
𝑋𝑏2 = ⎜0.2500
0.5625
⎝1.0000
0.0625⎟
0.3164
1.0000⎠
= � 0 2.2222 0 �
2.2222 0 2.0247
The variance-covariance matrix is
and the coefficient of biasing terms is
−3.8248
0.6404 0 −0.7031

𝑐𝑏 = �

�
5.0625
𝑀−1 (𝜉0 ) = �
0 0.4500 0 � ;
−0.7031 0 1.2656
The matrix of mean square error of each of the kth
segment is given by
𝑑𝑒𝑡𝑀(𝜉0 ) = 7.0233

−1 −1 ′

′ ′ −1

The standardized general variance function is given by
𝑀 = 𝑋

𝑘

+ 𝑋

𝑘

𝑋𝑘 𝑋𝑏𝑘 𝑐𝑏 𝑐𝑏 𝑋 𝑏𝑘 𝑋𝑘 𝑋

4

𝑔�𝑥� = 𝑛𝑥 ′ 𝑀−1 (𝜉0 )𝑥 = 𝛽0 + � 𝛽0 𝑥𝑖

2.7070 3.8821
 𝑀1 = �3.8821 6. � ; and

𝑖=1

𝑀 = � 2.7070
−3.8821
−3.8821
6.1347 �

IJSER © 2015 http://www.ijser.org

International Journal of Scientific & Engineering Research, Volume 6, Issue 1, January-2015 1851

ISSN 2229-5518

The matrices of convex combination of the segments
are 𝐻1 = 𝑑𝑖𝑎𝑔(0.5, 0.5) and 𝐻2 = (1 − 𝐻1 ) =
𝑑𝑖𝑎𝑔(0.5, 0.5). These matrices are normalized to have
𝐻∗ = 𝑑𝑖𝑎𝑔(0.7071, 0.7071), and
𝐻∗ = 𝑑𝑖𝑎𝑔(0.7071, 0.7071).
The direction vector 𝑑̂ , is 𝑑̂ = 𝑀𝑑 𝑧

𝑥1 = 0 = 𝑥𝑚𝑎𝑥.

Beginning another iteration in this iteration cycle, we
add the point, 𝑥𝑚𝑎𝑥 = 0, to the second segment (S2 ) and
repeat the process. The points in the first segment are
the same as to the first iteration, whereas the second segment is increased by the new point.
This is because; this point falls in the second segment.

𝑑 = �

𝑑0
� ; 𝑑̂ = �
0.6406
�
Thus, a MATLAB program was developed for this
𝑑1
27.5654
work and the result of the iteration cycles is present in
the table 5.1 below.


and the normalized direction vector, 𝑑̂ , such that 𝑑̂ ′

𝑑̂ = 1 is

𝑑̂ = �0.0232

0.9997
Then, the optimal starting point is

Table 5.1: The iteration sequence of D- optimal design in (-1, 1), for a quadratic single factor polynomial function

8

𝑥 = � 𝑤𝑖 𝑥𝑖

𝑖=1

1.0000

𝑥 = �

0

�
0.0000
In this table, the value of 𝜖 = 0.0005 is considered very
small, meaning that the optimal design is achieved.
Compute the optimal step-length, 𝜌0 , by substituting the
Therefore, the optimal design points are, (-1, 0.0048, 0,

values of 𝑥

0

= 0 and 𝑑̂ = 0.9997 in the variance
1).
function, we have
𝑔(𝜌0 ) = 2.5624 − 3.8248(0 + .9997𝜌0 )2 +
5.0624(0 + 0.9997𝜌0 )4

Problem 2

Let us consider the most general second-order polynomial in two factors with a known solution. Thus,

2 + 𝛽

𝑥2 ,

and solve for 𝜌0 to have
𝑓�𝑥� = 𝛽0 + 𝛽1 𝑥1 + 𝛽2 𝑥2 + 𝛽12 𝑥1 𝑥2 + 𝛽11 𝑥1

22 2

𝜌�0 = 0; 𝑜𝑟
𝜌�0 = 0.6148; 𝑜𝑟
𝜌�0 = 0.6148
−1 ≤ 𝑥1 , 𝑥2 ≤ 1; to obtain a 6-point exact D-optimal
design.
In the experimental region, 𝑋�, pick the following points.
Make a move to obtain a point,
�𝑥1
𝑥2
−1 −1 1
1 −1 −1
1 0
1 −1
−.25
�
. 25

𝑥1 = 𝑥0 + 𝜌�0𝑑̂

Then, the extended design matrix is
Thus, substituting the three values of 𝜌�0, and have

1 −1 1 −1 1 1

1 −1 −1 1 1 1

0 ⎛1 1

−1 −1 1 1 ⎞

𝑥1 = 0; 𝑜𝑟

𝑋(𝜉 ) = ⎜1 1 1 1

1 0 −1 0

1 1 ⎟

0 1

-0.6148; or

⎝1 −.25

. 25

−.0625

. 0625

. 0625⎠

0.6148.
Checking the prediction variance at each of these three points, we have
and the information matrix is

IJSER © 2015 http://www.ijser.org

International Journal of Scientific & Engineering Research, Volume 6, Issue 1, January-2015 1852

ISSN 2229-5518

𝑋 ′𝑋 = 𝑀(𝜉0 )

6 −.25 −.75

−.0625

4.0625

5.0025

1 −1 1
1 −1 0
1 1 1
1 1 0

⎛ −.25 4.0625 −0.625

. 0156

−.0156

−.0156

⎛1 −1
−1⎞ ; 𝑋 = ⎛1
1 −1⎞

= −.75

−.0625

−.0625

. 0156

5.0625

−.0156

−.0156

4.0039

. 0156

−.0039

−.9844⎞

−.0039⎟

𝑋1 = ⎜1 0
1 0
1 ⎟ 2
0
⎜1 0
1 0
−1⎟
0

4.0625

⎝ 5.0625

−.0156

−.0156

. 0156

−.9844

−.0039

−.0039

4.0039

4.00039

4.0039

5.0039⎠

⎝1 0
−1⎠
⎝1 0 1 ⎠
The variance-covariance matrix is
and their information matrices are respectively given as

𝑀−1 (𝜉0 )

1.1756 . 0667 −.0667

. 0167

. 05

−1.2422

𝑋�
6 −3 0

�2

6 3 0

⎛ . 0667 . 2500 0.0

0.0

0.0

−.0667 ⎞

1 = �−3 3 0� ; 𝑋

0 0 4
= �3 3 0� ;
0 0 4

= −.0667

. 0167

0.0

0.0

. 25

0.0

0.0

. 2500

−.25

0.0

. 3167

−.0167 ⎟ ; Where 𝑋�

1 1 2 2 2

. 0500

0.0

−.25

0.0

1.5000

−1.3

1 = 𝑋′ 𝑋

and 𝑋�
= 𝑋𝜄 𝑋

⎝−1.2422

−.0667

. 3167

−.0167

−1.3

2.5589 ⎠ The matrices of biasing terms, 𝑥1 𝑥2 , 𝑥1 , . . . , 𝑥2 , of the
𝐷𝑒𝑡𝑀 (𝜉0 ) = 225.000

2 4

variance function of each segment are
The standardized general variance function is

−1 1 1 −1 1

0 1 0 0 0

1 1 −1 −1 1 1

0 0 0 0 1 0

𝑑�𝑥, 𝜉0 � = 𝑛𝑥′ 𝑀−1 (𝜉0 )𝑥, = 𝑔(𝑥)

𝑋𝑏 = ⎛ 1

1 ⎜ 0

1 1 −1 −1 −1 1 1

0 1 0 0 1 0 0

1 1 1⎞

1 0 1⎟

𝑚 𝑚−1 𝑚 𝑚

2

= 𝛽0 + � 𝛽𝑖 𝑥𝑖 + � � 𝛽𝑖𝑗 𝑥𝑖 𝑥𝑗 + � 𝛽𝑖𝑖 𝑥𝑖
and

𝑖=1

𝑖=1 𝑖<𝑗

+. . .

𝑖=1

1 1 1 1 1

0 1 0 0 0

⎛

1 1 1 1 1 1

0 0 0 0 1 0

⎞

Where

′

= (1 𝑥 𝑥
𝑥 𝑥
𝑥2 𝑥2 )

𝑋𝑏2 = ⎜−1

1 1 1 −1

0 1 0 0

−1 1

−1 0

−1 −1 1 1

0 1 0 1⎟

𝑥1×𝑝

1 2 1 2 1 2 ⎠

Thus, the variance function is
𝑔�𝑥� = 7.0536 + .80042𝑥1 − .8004𝑥2 + .2004𝑥1 𝑥2
while the coefficient of biasing terms is

𝑐′ = (.2004 2.1000 − 13.4064 − .8004 − 3.000 −

2 2 2

+ 2.1𝑥1 − 13.4064𝑥2 − .8002𝑥1 𝑥2

3.8004 − 14.1000 − .002 − .1002 9.0000 15.3534 )

2 3 2 2

− 3𝑥1 𝑥2 + 3.8004𝑥2 − 14.1𝑥1 𝑥2
− .002𝑥1 𝑥3 − .1002𝑥3 𝑥 + 9𝑥4

2 1 2 1

+ 15.3534𝑥4
The matrix of mean square error of each of the kth
segment is given by

−1 −1 ′

′ ′ −1

We now partition the experimental region −1 ≤
𝑥1 , 𝑥2 ≤ 1 into k-segments, k=2 with common
𝑀 = 𝑋

𝑘

+ 𝑋

𝑘

𝑋𝑘 𝑋𝑏𝑘 𝑐𝑏 𝑐𝑏 𝑋 𝑏𝑘 𝑋𝑘 𝑋

boundary as
𝑆1 = {−1 ≤ 𝑥1 ≤ 0, −1 ≤ 𝑥2 ≤ 1} and 𝑆2 =
{0 ≤ 𝑥1 ≤ 1, −1 ≤ 𝑥2 ≤ 1}, and in the kth segment,
𝑔�𝑥� is approximated by a first-order linear function.
Thus, 𝑦𝑘 = 𝛽00 + 𝛽1𝑘 𝑥1+ 𝛽2𝑘 𝑥2; k=1, 2.
Pick points in the 𝑆𝑘 segments, and obtain the 𝑋𝑘
design matrices as thus

2.0181 −2.5659 2.9859

=> 𝑀1 = �−2.5659 5.6556 −5.1382� ; 𝑎𝑛𝑑

2.9895 −5.1382 5.5418

2.0181 1.1807 2.9895

The matrices of convex combination of the segments are arranged in the Hi matrices as follows
. 5000 0 0
𝐻1 = �
0 . 2639 0 � ;
0 0 . 5000

IJSER © 2015 http://www.ijser.org

International Journal of Scientific & Engineering Research, Volume 6, Issue 1, January-2015 1853

ISSN 2229-5518

0 0 . 5000
The point of maximum prediction variance among the three points is
and the normalized Hi are:
𝑥∗ = �
−.0180
�
−.0344
. 7071 0 0

∗

Add this point 𝑥1 , to the second segment (𝑆2 ) of the
𝐻1 = �
0 . 3374 0 � ;
experimental region and repeat the process. Notice that
0 0 . 7071
. 7071 0 0
the points in the first segment have increased by the
addition of the new point, whereas the points in the
𝐻2 = �
0 . 9414 0 �
second segment remain the same.
0 0 . 7071
The direction vector 𝑑̂ , is

𝑑̂ = 𝑀−1 𝑧

𝑑

−0.1069
= � 9.4973 �
The summary of the iteration cycles is given in the table below.

Table 5.2: The iteration sequence of D-optimal design in [-1, 1], for a quadratic polynomial function of two factors

−52.6626

and normalized the direction vector, 𝑑̂ ,
such that 𝑑̂ ′ 𝑑̂ = 1 is

𝑑̂ = �−.0020�

. 1775
−.9841
Then, the optimal starting point is

12

𝑥 = � 𝑤𝑖 𝑥𝑖

𝑖=1

1.0000

𝑥 = � �

0 −0.0031

−0.1166
The value of 𝜖 = .0003 is considered very small, meaning that the optimal design is reached.
Atkinson and Donev (1992) used KL exchange method to solve the problem for generating an exact D-optimal design for this regression model and obtain the same result which Box and Draper (1971) obtained via a computer hill-climbing search.
Their D-optimal design is

𝑛

= (−1, −1), (1, −1), (−1,1), (−.1315, −.1315), (. 3945, 1), (1, .3945)

Compute the optimal step-length 𝜌0, by substituting the

values of 𝑥

0

and 𝑑̂ in the variance function, and

5.2 Assessment of the search method

thereafter differentiate with respect to 𝜌0, equate to zero and solve for 𝜌�0, to obtain
𝜌0 = .6483, 𝑜𝑟 𝜌�0 = −.7296; 𝑜𝑟 𝜌�0 = −.0835
Make a move to obtain a point

𝑥1 = 𝑥0 + 𝜌�0𝑑̂

Then, substituting the value of 𝜌�0, we have
The assessment of this method was considered in two folds. The first one is the computer execution
time (c. e. t.) of the algorithm to get to D-optimal design. In this paper, a Matlab program of this algorithm was developed and average of twenty running time of each model is presented in the table below. All runs were performed in a Presario CQ56 laptop computer.
. 1119
� ; 𝑜𝑟 𝑥
−.1326
= � � ; 𝑜𝑟 𝑥1
The second other way of assessment is through D- efficiency. D-efficiency is a relative measure of how

𝑥1 = �−.7546

1 . 6015

−.0180
= � �
−.0344
design compares with the D-optimum design for a
specific design experimental region. Design efficiency

IJSER © 2015

http://www.ijser.org

International Journal of Scientific & Engineering Research, Volume 6, Issue 1, January-2015 1854

ISSN 2229-5518

is used to compare user-specified design to the optimal design. Theoretically, design efficiency lies between 0 and 1, and the closer the efficiency to 1, the better the arbitrary design. If information matrix for optimal
design is 𝑀(𝜉1 ), and suppose an arbitrary information matrix of a design is 𝑀(𝜉2 ). Then the D-efficiency of
the arbitrary design is defined as

|𝑀(𝜉2 )| 1

method is a single search system, in contrast to the “Adjustment Algorithm” of Atkinson and Donev (1992, Chapter 15), where a secondary search is required at the end of the primary one before the search terminates.

References

[1] Atkinson, A. C. and Donev, A.C. (1989). The
Construction of D-optimal Experimental Exact

𝐷𝑒𝑓𝑓 = {
}𝑝 , where p is the number of model
Designs with Application to Blocking

|𝑀(𝜉1 )|

parameters: see Atkinson and Donev (1989).
The table below shows how the method compares with the KL-exchange method reported by Atkinson and Donev.

Table 5.3: The performance of method

The table 5.3 above shows that the method is efficiently equal to what Atkinson, Donev and Tobias reported for a single factor quadratic polynomial function. In the same way, the new method is highly efficient as the result reported using KL exchange method for two factor quadratic polynomial models. Again, the new method reached the optimal design points in only four iteration cycles in each of the example with minimal execution time.

6. Conclusion

A line search method for constructing D-optimal

exact designs for a polynomial response function f(𝑥)

of degree z, where the factor levels 𝑥 are defined in
continuous geometry is shown. Under this search
method, variance of the regression function is approximated by a polynomial function of degree 2z. A maximizer of this polynomial is determined by segmented line search technique. Thereafter, this maximizer is used to replace the point of minimum prediction variance from amongst the point already in the design. This exchange of points is continued until the sequence converges to D-optimal design.
The search method compares favorably with the method of KL-Exchange Algorithm reported by Atkinson and Donev. We affirm that the new search
Response Surface Designs. Biometrics, 176,
515-527
[2] Atkinson, A. C and Donev, A.C. (1992). Optimum Experimental Designs. New York, Oxford University Press.
[3] Atkinson, A. C. Donev, A.C., and Tobias, R.D. (2006). Optimum Experimental Designs, with SAS. London, Macclesfield and Cary
[4] Chigbu, P.E and Oladugba, A.V. (2010). On the Loss of Information Matrix and Its Analytical Properties in DOE (Design of Experiments). J. of Nigerian Math. Society. Vol. 29, 41-49
[5] Fedorov, V. V. (1972). Theory of Optimal
Experiments. Academic Press, New York.
[6] Feller, W. (1966). An Introduction to Probability Theory and its Applications, 3rd Edition, John Wiley and Sons Inc., New York.
[7] Gaffke, N. and Krafft, O. (1982). Exact D-Optimal Designs for Quadratic Regression. Journal of the Royal Statistical Society, Series B V(44), N(2). 394-397.
[8] Onukogu, I.B. (1997). Foundation of Optimal Exploration of Response Surface. Ephrata Nsukka.
[9] Onukogu, I.B. and Chigbu, P.E. (2002). Super Convergent Line Series in Optimal Design of Experiment and Mathematical Programming, AP Express Publishers, Nsukka. K
[10] Onukogu, I.B. and Nsude F.I. (2014), A global convergent sequence for construction of D- optimal exact designs in an infinite space of trials. Journal of Statistics: Advances and Theory of Applications, 12(1), pp 39-52
[11] Ugbe, T.A. and Chigbu, P.E. (2014). On Non Overlapping Segmentation of the Response Surface for Solving Constrained Programming Problems through Super Convergent line

IJSER © 2015 http://www.ijser.org

International Journal of Scientific & Engineering Research, Volume 6, Issue 1, January-2015 1855

ISSN 2229-5518

Series. Communication in Statistics-Theory and Methods,Vol.43,306-320.

IJSER © 2015 http://www.ijser.org