# Fast Computation of Highly $G$-Optimal Exact Designs via Particle Swarm
Optimization

Stephen J. Walsh and John J. Borkowski

###### Abstract

Computing proposed exact $G$-optimal designs for response surface models is a
difficult computation that has received incremental improvements via algorithm
development in the last two-decades. These optimal designs have not been
considered widely in applications in part due to the difficulty and cost
involved with computing them. Three primary algorithms for constructing exact
$G$-optimal designs are presented in the literature: the coordinate exchange
(CEXCH), a genetic algorithm (GA), and the relatively new $G$-optimal via
$I_{\lambda}$-optimality algorithm ($G(I_{\lambda})$-CEXCH) which was
developed in part to address large computational cost. Particle swarm
optimization (PSO) has achieved widespread use in many applications, but to
date, its broad-scale success notwithstanding, has seen relatively few
applications in optimal design problems. In this paper we develop an extension
of PSO to adapt it to the optimal design problem. We then employ PSO to
generate optimal designs for several scenarios covering $K=1,2,3,4,5$ design
factors, which are common experimental sizes in industrial experiments. We
compare these results to all $G$-optimal designs published in last two decades
of literature. Published $G$-optimal designs generated by GA for $K=1,2,3$
factors have stood unchallenged for 14 years. We demonstrate that PSO has
found improved $G$-optimal designs for these scenarios, and it does this with
comparable computational cost to the state-of-the-art algorithm
$G(I_{\lambda})$-CEXCH. Further, we show that PSO is able to produce equal or
better $G$-optimal designs for $K=4,5$ factors than those currently known.
These results suggest that PSO is superior to existing approaches for
efficiently generating highly $G$-optimal designs.

##  1 Introduction

The groundbreaking publication which defined the optimal design concept
predates R. A. Fisher’s foundational 1935 work _The Design of Experiments_
(atkRev; fisher:1935; firstGopt). Kirstine Smith, a Danish female statistician
working under Karl Pearson, published guidance in a 1918 Biometrika paper on
how the residual error variance suggests designing an experiment used to fit a
polynomial (firstGopt). atkRev note that this contribution is considered the
seminal optimal design paper and was 30 years ahead of its time. This paper’s
contribution provided what we now refer to as $G$-optimal designs for a $K=1$
factor experiment supporting a first up to sixth-degree polynomial.

The exact $G$-optimal design generation problem is recognized as the
notoriously difficult mini-max problem: we are seeking designs that minimize
the maximum prediction variance over the study space. It wasn’t until recent
decades that modern computing resources were cheap enough to enable the
application of well-suited meta-heuristic optimization approaches to this
problem. jobo1 provided one of the earliest applications of Genetic Algorithms
(GAs) for generating _exact_ optimal designs over the hypercube supporting a
second-order response surface model in $K=1,2,3$ design factors and for
experiment sizes $N=3,4,5,6,7,8,9$, $N=6,7,8,9,10,11,12$, and
$N=10,11,12,13,14,15,16$ respectively. jobo1 adapted a GA to generate $G$-,
$D$-, $I$-, and $A$-optimal designs and published a proposed catalog of exact
optimal designs for each criterion. In contrast to continuous designs which
can appeal to the _General Equivalence Theorem_ to verify the global
optimality of a proposed continuous design, no such general theory exists for
exact designs and so the community recognizes that proposed global exact
optimal designs must be published and then challenged by future developments
and algorithm applications. Therefore, the results provided in jobo1 have
become a ‘ground truth’ standard data set against which to compare results of
newly developed algorithms to solve this problem. Until this point, the exact
$G$-optimal designs of Borkowski (2003) remained the best-known designs for
these scenarios and were subsequently reproduced by rodman via an augmented
application of the coordinate exchange and used as a benchmark by gIopt.

The next large investigation into generating and establishing the properties
of exact $G$-optimal designs is provided by rodman. In this paper, the authors
adapted the coordinate exchange (CEXCH) algorithm of mncoord in conjunction
with Brent’s minimization algorithm (used to score a candidate design on the
$G$-scale) to generate exact $G$-optimal designs. These authors covered all
scenarios in Borkowski (2003) and were able to reproduce these designs with
their new algorithm. A second contribution of this work was an extension of
the algorithm to generate designs for a few $K=4,5$ factor cases—a design
scenario not addressed to that point due to computation cost. The authors
propose candidate exact $G$-optimal designs for $K=4$, $N=15,20,24$ and $K=5$,
$N=21,26,30$ which support the second-order response surface model. A third
contribution of this work was a detailed comparison of the resulting $G$-CEXCH
generated exact $G$-optimal design’s properties (via fraction of design space
plots) with corresponding $D$\- and $I$-optimal designs. The authors found
that, for scenarios where post-hoc prediction was the primary objective, the
$I$-optimal designs exhibited smaller prediction variance than the
corresponding $G$-optimal designs over large percentages of the design space
(even though the $G$-optimal design has lower maximum prediction variance).
This observation motivated the authors to recommend that the additional
computing cost required for $G$-optimal designs may not be worth the effort
and $I$-optimal designs might be preferred for many practical experimental
scenarios. Nonetheless, this catalog of proposed exact $G$-optimal designs for
small response surface models remains a useful and referenced data set against
which computational statisticians can compare the performance of the
adaptation of new algorithms to this problem.

The most recent publication on generating exact $G$-optimal designs is offered
by gIopt. The primary focus of this paper is the computational cost associated
with generating $G$-optimal designs. These authors noticed a correspondence
between continuous $I_{\lambda}$-optimal designs and $G$-optimal designs. They
propose a new algorithm that exploits the structural relationship between
these two types of designs which first generates a continuous
$I_{\lambda}$-optimal design as a starting point for the exact $G$-optimal
design search, Then they apply the CEXCH algorithm to locally improve the
candidate $G$-optimal design. We will refer to their algorithm as
$G(I_{\lambda})$-CEXCH. These authors are the first to completely run
$G(I_{\lambda})$-CEXCH and rerun the searches via GA by Borkowski (2003) in
Matlab, and CEXCH by rodman in JMP (they also reported the architecture of
their PC) in order to compare computing efficiency of their proposed algorithm
to the others. They implemented 200 runs of each algorithm for each design
scenario discussed in jobo1 and also ran two new scenarios: $K=4$, $N=17$ and
$K=5$, $N=23$ with only $G(I_{\lambda})$-CEXCH (they indicated that $G$-CEXCH
would have taken 25 and 166 days respectively to complete the task, and so
this algorithm was not run on these scenarios). To briefly summarize their
results, they found that, with respect to finding the highest $G$-optimal
designs, the GA was superior (finding the best design each set of 200
algorithm runs) but that $G$-CEXCH and $G(I_{\lambda})$-CEXCH produced designs
with high relative efficiency vs. $GA$ generated designs with relative-
efficiencies of 90% or greater. However, the GA achieves this superiority at
high computational cost exhibiting 2-orders of magnitude increase in the
number of objective function evaluations as compared to the
$G(I_{\lambda})$-CEXCH. $G(I_{\lambda})$-CEXCH was demonstrated to produce
designs with high $G$-optimality as well as $G$-CEXCH but with computing times
at a fraction of $G$-CEXCH for the higher dimensional searches. Thus they
propose that $G(I_{\lambda})$-CEXCH is a good choice for generating exact
$G$-optimal designs as it is 1.) demonstrated to produce designs with high
relative efficiency to GA generated designs and 2.) does this more efficiently
than existing algorithms.

Particle swarm optimization (PSO) is a mature, simple, but highly effective
meta-heuristic optimization algorithm (well suited for finding global optima
on multimodal objectives) in the class of evolutionary algorithms (including
GA) and has not been applied to this problem to date. Since its publishing by
opso_1995 in the machine learning literature, this algorithm has attained
widespread application and success in many scientific domains—to date Google
tracks this paper with over 69,216 citations (queried January 17, 2022) which
illustrates its widespread impact in algorithm development research and
applications to engineering and applied science domains. This success
notwithstanding, PSO has seen relatively few applications in optimal design,
and for those publications that do address the topic, they tend to focus on
using PSO to generate optimal _continuous_ designs—see, for example, chen0;
wong; chen1; chen2; Chen3; chen4. There have also been some applications of
PSO, with apparent high success, outside of the response surface design
framework—see, for example, lukemire for an application of PSO to generating
$D$-optimal designs for logistic-regression models and joseph for an
application of PSO to generating space filling designs. The contribution of
our paper is: 1) a first application of PSO to the exact $G$-optimal design
generation problem, 2) a review and summary of PSO generated results against
those algorithms and results published over the last 20 years, and 3) we
update the proposed catalog of exact $G$-optimal designs with new-found
improved $G$-optimal designs for many of the studied scenarios.

In the next section we present our notation and define the exact $G$-optimal
design generation problem. In section 3 we develop a yet unpublished extension
of PSO for generating exact optimal designs. In section 4 we compare $G$-PSO
generated designs to those provided by jobo1, rodman and gIopt. In section 5
we provide discussion, conclusions, and future research directions. Based on
the new information reported in this paper, we conclude with a recommendation
for which algorithm appears to be best suited for generating exact $G$-optimal
designs.

##  2 $G$-Optimal Designs

As in previous publications, we will be working with the second-order response
surface model under standard assumptions. Our notation is as follows. Let
$\mathcal{X}$ denote a design space. Without loss of generality with respect
to generating optimal designs on irregular hyper-rectangular design regions,
it is common, in a response surface setting, to assume that the design factors
have been scaled to $[-1,1]$. In which case the design space is the
$K$-dimensional hypercube $\mathcal{X}=[-1,1]^{K}$ for $K$ factors. Let
$\mathbf{x}^{\prime}\in\mathcal{X}$ denote a $1\times K$ point in the design
space. For $N$ specified design points, the set of distinct design points
$\mathbf{x}^{\prime}_{i},\text{ }i=1,2,\ldots,N$ are organized as the rows of
the $N\times K$ design matrix

$\mathbf{X}=\begin{pmatrix}\mathbf{x}^{\prime}_{1}\\\
\mathbf{x}^{\prime}_{2}\\\ \vdots\\\ \mathbf{x}^{\prime}_{N}\end{pmatrix}.$

Next, we provide some notation that clarifies the difference between the
design space and the model space—the latter of which is required to evaluate
the properties and the quality of a candidate design $\mathbf{X}$. Let $p$
represent the number of parameters of the model that one intends to fit to the
data. We denote the $N\times p$ model matrix as $\mathbf{F}(\mathbf{X})$. Here
we have emphasized that the model matrix is a function of the design matrix.
For simplicity, however, from this point forward we will simply represent the
model matrix by the shorthand $\mathbf{F}$. A similar device is required for a
model vector, which is a function that maps a candidate design point to the
model space. We adopt the notation $\mathbf{f}^{\prime}(\mathbf{x}^{\prime})$
which is a $1\times p$ row-vector.

We provide an example for the case of an $N$ point design with $K=3$ design
factors on the second order model to clarify the use of this notation. In this
case there are

$p={K+2\choose 2}={5\choose 2}=10$

parameters in the second-order model (rsm; jobo1). The design matrix for this
model is

$\mathbf{X}=\begin{pmatrix}\mathbf{x}^{\prime}_{1}\\\
\mathbf{x}^{\prime}_{2}\\\ \vdots\\\
\mathbf{x}^{\prime}_{N}\end{pmatrix}=\begin{pmatrix}x_{11}&x_{12}&x_{13}\\\
x_{21}&x_{22}&x_{23}\\\ &\vdots&\\\ x_{N1}&x_{N2}&x_{N3}\end{pmatrix}.$

The second-order model, for a single design point
$\mathbf{x}^{\prime}=(x_{1}\quad x_{2}\quad x_{3})$ in scalar form is written

$y=\beta_{0}+\sum_{i=1}^{K}\beta_{i}x_{i}+\sum_{i=1}^{K-1}\sum_{j=i+1}^{K}\beta_{ij}x_{i}x_{j}+\sum_{i=1}^{K}\beta_{ii}x_{i}^{2}+\epsilon,$

and so the mapping of the design point $\mathbf{x}^{\prime}$ into the model
space is the $1\times p$ vector

$\mathbf{f}^{\prime}(\mathbf{x}^{\prime})=\begin{pmatrix}1&x_{1}&x_{2}&x_{3}&x_{1}x_{2}&x_{1}x_{3}&x_{2}x_{3}&x_{1}^{2}&x_{2}^{2}&x_{3}^{2}\end{pmatrix}.$

The mapping of the design matrix into the model matrix is

$\mathbf{F}(\mathbf{X}):=\mathbf{F}=\begin{pmatrix}1&x_{11}&x_{12}&x_{13}&x_{11}x_{12}&x_{11}x_{13}&x_{12}x_{13}&x_{11}^{2}&x_{12}^{2}&x_{13}^{2}\\\
1&x_{21}&x_{22}&x_{23}&x_{21}x_{22}&x_{21}x_{23}&x_{22}x_{23}&x_{21}^{2}&x_{22}^{2}&x_{23}^{2}\\\
\vdots&&\vdots&&\vdots&&\vdots&&\vdots&\\\
1&x_{N1}&x_{N2}&x_{N3}&x_{N1}x_{N2}&x_{N1}x_{N3}&x_{N2}x_{N3}&x_{N1}^{2}&x_{N2}^{2}&x_{N3}^{2}\end{pmatrix}.$

This enables writing the full second-order linear model in the $N$ design
points in matrix form

$\mathbf{y}=\mathbf{F}\bm{\beta}+\bm{\epsilon},$

where $\bm{\beta}$ is a $p\times 1$ vector of parameters. Under the ordinary
least squares (OLS) framework, it is assumed that
$\bm{\epsilon}\sim\mathcal{N}_{N}(\mathbf{0},\sigma^{2}\mathbf{I}_{N})$ where
$\mathcal{N}_{N}$ denotes the $N$-dimensional multivariate normal
distribution. The least-squares estimate of the $\bm{\beta}$’s (which in this
framework are also ML estimates) is

$\hat{\bm{\beta}}=(\mathbf{F}^{\prime}\mathbf{F})^{-1}\mathbf{F}^{\prime}\mathbf{y},$

which have variance-covariance

$\text{Var}(\hat{\bm{\beta}})=\sigma^{2}(\mathbf{F}^{\prime}\mathbf{F})^{-1}.$

A prediction of the response at new point $\mathbf{x}^{\prime}$ is denoted as
$\hat{y}(\mathbf{x}^{\prime})$, and the variance of the mean predicted value
is

$\text{Var}(\hat{y}(\mathbf{x}^{\prime}))=\sigma^{2}\mathbf{f}^{\prime}(\mathbf{x}^{\prime})(\mathbf{F}^{\prime}\mathbf{F})^{-1}\mathbf{f}(\mathbf{x}^{\prime}).$

The total information matrix for $\bm{\beta}$, specifically
$\mathbf{M}(\mathbf{X})=\mathbf{F}^{\prime}\mathbf{F}$, plays an important
role in optimal design of experiments—all optimal design objective functions
are functions of this matrix.

While $\mathcal{X}$ denotes the space of candidate _design points_
$\mathbf{x}^{\prime}$, a _design matrix_ $\mathbf{X}$ is a collection of $N$
such design points. Thus, the space of all candidate designs can be denoted

$\mathbf{X}\in\bigtimes_{j=1}^{N}\mathcal{X}=\bigtimes_{j=1}^{N}[-1,1]^{K}=[-1,1]^{NK}=\mathcal{X}^{N},$
(1)

and the practitioner must choose a design from $\mathcal{X}^{N}$ to implement
the experiment in practice. Note that the optimization problem has dimension
$N*K$ (the dimension of the optimization search for the designs studied here
is sometimes, erroneously, referred to as ‘small’ and only viewed as a
function of the number of design factors $K$). An optimality criterion is used
to define which candidate designs $\mathbf{X}\in\mathcal{X}^{N}$ are ‘good’
designs. An optimization algorithm (such as GA, $G$-CEXCH, or
$G(I_{\lambda})$-CEXCH) is used to search $\mathcal{X}^{N}$ to find the
‘best’, or hopefully, globally optimal design. A $G$-optimal design is that
design which minimizes the maximum scaled prediction variance over all points
of prediction in $\mathcal{X}$. Recall that the variance of a mean prediction
is

$\text{Var}(\hat{y}(\mathbf{x}^{\prime}))=\sigma^{2}\mathbf{f}^{\prime}(\mathbf{x}^{\prime})(\mathbf{F}^{\prime}\mathbf{F})^{-1}\mathbf{f}(\mathbf{x}^{\prime})$

for any point of prediction $\mathbf{x}^{\prime}\in\mathcal{X}$. The scaled
prediction variance (SPV) removes the scale parameter $\sigma^{2}$ and re-
scales to $N$ by multiplying by the factor $N/\sigma^{2}$. Thus, SPV is
defined as (for a candidate design $\mathbf{X}$)

$\displaystyle\text{SPV}(\mathbf{x}^{\prime})$
$\displaystyle:=\mathbf{f}^{\prime}(\mathbf{x}^{\prime})\mathbf{M}^{-1}(\mathbf{X})\mathbf{f}(\mathbf{x}^{\prime})$
$\displaystyle=N\mathbf{f}^{\prime}(\mathbf{x}^{\prime})(\mathbf{F}^{\prime}\mathbf{F})^{-1}\mathbf{f}(\mathbf{x}^{\prime}).$
(2)

Note that defining $SPV$ as in the second line of Eq. 2 is necessary, as will
be explained leading up to Eq. 5.

The $G$-score of an arbitrary design $\mathbf{X}$ is defined as the maximum
scaled prediction variance over all points of prediction
$\mathbf{x}^{\prime}\in\mathcal{X}$

$\displaystyle G(\mathbf{X})$
$\displaystyle:=\max_{\mathbf{x}^{\prime}\in\mathcal{X}}\text{SPV}(\mathbf{x}^{\prime})$
$\displaystyle=\max_{\mathbf{x}^{\prime}\in\mathcal{X}}N\mathbf{f}^{\prime}(\mathbf{x}^{\prime})(\mathbf{F}^{\prime}\mathbf{F})^{-1}\mathbf{f}(\mathbf{x}^{\prime}).$
(3)

Eq. 2 indicates that for a fixed candidate design $\mathbf{X}$, scoring the
candidate on the $G$ scale is itself an optimization problem. The $G$-optimal
design $\mathbf{X}^{*}$ is that design which minimizes, over all designs
$\mathbf{X}\in\mathcal{X}^{N}$, the maximum scaled prediction variance, namely

$\displaystyle\mathbf{X}^{*}$
$\displaystyle:=\operatorname*{argmin}_{\mathbf{X}\in\mathcal{X}^{N}}G(\mathbf{X})$
$\displaystyle=\operatorname*{argmin}_{\mathbf{X}\in\mathcal{X}^{N}}\max_{\mathbf{x}^{\prime}\in\mathcal{X}}\text{SPV}(\mathbf{x}^{\prime})$
$\displaystyle=\operatorname*{argmin}_{\mathbf{X}\in\mathcal{X}^{N}}\max_{\mathbf{x}^{\prime}\in\mathcal{X}}N\mathbf{f}^{\prime}(\mathbf{x}^{\prime})(\mathbf{F}^{\prime}\mathbf{F})^{-1}\mathbf{f}(\mathbf{x}^{\prime}).$
(4)

Eq. 2 shows that finding the $G$-optimal design is a minimax problem. This
optimization has proved notoriously difficult to solve because neither of the
optimizations required to compute $\mathbf{X}^{*}$ are convex. The scale of
$G$ for an arbitrary design has known bounds. The General Equivalence Theorem
of kiefer demonstrates that the lower bound on $G(\mathbf{X})$ is

$G(\mathbf{X})=\max_{\mathbf{x}^{\prime}\in\mathcal{X}}\text{SPV}(\mathbf{x}^{\prime})\geq
p.$ (5)

That is, the smallest that the maximum scaled prediction variance of candidate
design $\mathbf{X}$ may be is $p$, the number of parameters. This apprently
gives a way to verify if a proposed exact $G$-optimal design is globally
optimal, however, not all design scenarios have globally $G$-optimal designs
that will achieve this lower bound. Further, if design $\mathbf{X}^{*}$ is
globally $G$-optimal and achieves the result in Eq. 5, then for this design
$\text{SPV}(\mathbf{x^{\prime}})=p$ at all diagonals of the hat matrix
$\mathbf{F}(\mathbf{F}^{\prime}\mathbf{F})^{-1}\mathbf{F}^{\prime}$ and
$\text{SPV}(\mathbf{x^{\prime}})\leq p$ at all other points of prediction
$\mathbf{x}^{\prime}\in\mathcal{X}$ (rsm). These are the only results from the
theory of continuous optimal designs that also apply directly to exact optimal
designs. The result stated in Eq. 5 is an important result—it means that the
scale of $G$-optimality is not arbitrary (at least up to the number of
parameters in the model $p$). It is customary to exploit this fact and score
candidate designs on the $G$-efficiency scale,

$G_{\text{eff}}(\mathbf{X})=100\times\frac{p}{G(\mathbf{X})},$ (6)

in order to gauge the quality of a candidate design $\mathbf{X}$ (larger
$G_{\text{eff}}$ on this scale implies a better design). Last, the _relative
efficiency_ of two candidate designs may be computed as

$G_{\text{releff}}(\mathbf{X}_{1},\mathbf{X}_{2})=100\frac{G_{\text{eff}}(\mathbf{X}_{1})}{G_{\text{eff}}(\mathbf{X}_{2})}$
(7)

The reason for constructing an algorithm to search $\mathcal{X}^{N}$ globally
for the $G$-optimal design is that once one has a candidate design
$\mathbf{X}$, one must gauge its quality via Eq. 2. In the context of the
larger search of $\mathcal{X}^{N}$ for the globally optimal design
$\mathbf{X}$, each candidate $\mathbf{X}$ must be scored on the $G$-scale via
Eq. 2, and this is a separate optimization of a multimodal function. There are
two published approaches in the literature suggesting ways to score a
candidate matrix $\mathbf{X}$ via computing or approximating $G(\mathbf{X})$.

First, authors jobo1 and gIopt exploit the symmetry of the SPV-surface for a
$G$-optimal design and recommend using a $5^{K}$ point grid over $\mathcal{X}$
with each factor containing $x_{j}\in\\{-1,0.5,0,0.5,1\\}$ for $j=1,\ldots,K$
grid-points. The full grid defined as
$\mathbf{G}_{\mathcal{X}}=\\{-1,0.5,0,0.5,1\\}^{K}$ (gIopt). Once one has a
candidate design, the SPV is evaluated on $\mathbf{G}_{\mathcal{X}}$ and the
maximum value is taken as an approximation of the $G$-score for the candidate
design. While this approach computes an approximate $G$-score for a candidate
design $\mathbf{X}$ both authors have noted that, due to the symmetry of the
$G$-surface for $G$-optimal designs, the approximation is quite adequate and
supports finding highly $G$-optimal designs.

A second approach is published by authors rodman who apply a different
approach to scoring a candidate design on the $G$-scale during the evolution
of the $G$-CEXCH algorithm. In this algorithm once a candidate is available
(e.g. after every coordinate exchange) they employ Brent’s minimization
algorithm to search for the maximum prediction variance for the candidate
design over

Regarding the choice of method for scoring a candidate design on the $G$-scale
in this paper: we implemented the $5^{K}$ grid approach recommended by jobo1
and gIopt for several reasons. First, it is easy to implement in the context
of the larger $G$-optimal design search on $\mathcal{X}^{N}$. Second, in
preliminary studies we also found these approximations to the $G$-score to be
adequate and subject to small error. Third, in preliminary studies, we
attempted a separate optimization approach similar to rodman but found that if
a candidate design was mis-scored the effect was disastrous on the evolution
of the algorithm searching $\mathcal{X}^{N}$ for the optimal design.

In the next section we describe the PSO algorithm and provide an extension of
PSO to generate $G$-optimal designs used to produce the results presented in
this paper.

##  3 Particle Swarm Optimization

PSO has been demonstrated to perform very well for high-dimensional
optimization of multimodal objective functions, see for example, opso_1995;
spso1_2007; spso2_2011; hal_spso; opso2_1995; ken_1997; eberhart_1998; clerc2;
engel; zombook; bratton. PSOs strengths include: 1) few-to-no assumptions
about the properties of the objective function $f$ to be optimized, 2) PSO is
demonstrated to be strongly robust to entrapment in local optima, and thereby
is a good ‘match’ to the exact optimal design generation problem, 3)
simplicity—the core function of the algorithm can be explained via two simple
update equations, and 4) in contrast to other meta-heuristics where studying a
range of tuning parameters can yield more efficient searches for specific
problems, PSO only has three tuning parameters and these have been studied
extensively theoretically and empirically with optimal values demonstrated for
searches (such as ours) that fall on the typical Euclidean geometry, see
clerc2; ebershi_3; clerc1; shi2 among others. Therefore, we employ the values
of the PSO tuning parameters as suggested in the machine learning literature,
and these values are guaranteed to 1) ensure a large degree of swarm
exploration of the search space and 2) ensure that the particles in the swarm
eventually settle on a consensus point for the optimal solution to the
optimization problem (clerc2; hal_spso).

PSO is suited to perform the calculation

$\mathbf{x}^{*}=\operatorname*{argmin}_{\mathbf{x}\in\mathcal{X}}f(\mathbf{x})$
(8)

for objective $f$ with $K\times 1$ input vector $\mathbf{x}$. The algorithm
proceeds by first randomly drawing several candidate solution points (i.e. the
particles) $\mathbf{x}_{i}$ according to a multivariate uniform distribution
in $\mathcal{X}$ for $i=1,\ldots,S$ particles. Similarly each particle is
assigned an initial velocity $\mathbf{v}_{1}$ which governs particle $i$’s
step size and direction for its search of the optimizer of $f$. Then, for
iteration $t$ of the PSO algorithm the particles step through the search space
(searching for the optimum of $f$) according to velocity and position update
equations

$\displaystyle\mathbf{v}_{i}(t+1)$ $\displaystyle=\omega\mathbf{v_{i}}(t)$
(inertia) (9) $\displaystyle+$ $\displaystyle
c_{1}U_{K}(\mathbf{0},\mathbf{1})\odot\left(\mathbf{p}_{\text{best},i}-\mathbf{x}_{i}(t)\right)$
(cognitive) (10) $\displaystyle+$ $\displaystyle
c_{2}U_{K}(\mathbf{0},\mathbf{1})\odot\left(\mathbf{g}_{\text{best}}-\mathbf{x}_{i}(t)\right)$
(social) (11) $\displaystyle\mathbf{x}_{i}(t+1)$
$\displaystyle=\mathbf{x}_{i}(t)+\mathbf{v}_{i}(t+1)$ (position update) (12)

where $\odot$ represents the Hadamard product (elementwise multiplication),
$\mathbf{p}_{\text{best},i}$ is the ‘personal best’ position in $\mathcal{X}$
regarding fitness on $f$ that particle $i$ has visited to time $t$ (e.g. this
is the particle’s ‘memory’), $\mathbf{g}_{\text{best}}$ represents the ‘group
best’ position in $\mathcal{X}$ regarding fitness on $f$ that particle $i$
found by the swarm to time $t$ (e.g. this the the swarm’s ‘communication’),
and $U_{K}(\mathbf{0},\mathbf{1})=\\{u_{j}\\}$ for $j=1,\ldots,K$ where
$u_{j}\overset{i.i.d.}{\sim}U(0,1)$ is a random scaling factor for each
velocity update component. Finding optimal velocity weighting constants
$\omega,c_{1},$ and $c_{2}$ was the subject of intense research over
approximately a decade. In this work we employ, as recommended by the
extensive PSO literature for optimization problems of our nature,
$\omega=0.72984,c_{1}=c_{2}=2.05\omega=1.496172$. These values are guaranteed
to balance particle exploration while ensuring that the velocities eventually
scale down so that the swarm radius eventually decreases to 0, i.e. the
particles in the swarm come to a consensus regarding a solution to
$\operatorname*{argmin}_{\mathbf{x}\in\mathcal{X}}f(\mathbf{x})$ (hal_spso;
clerc1; clerc2). Last, the personal best position of particle $i$,
$\mathbf{p}_{\text{best},i}$, and the group’s best position,
$\mathbf{g}_{\text{best}}$, are updated at each iteration $t$ via the simple
logic

$\displaystyle\textbf{if}f(\mathbf{x}_{i})<f(\mathbf{p}_{\text{best},i})$
$\displaystyle\hskip
11.38092pt\mathbf{p}_{\text{best},i}\leftarrow\mathbf{x}_{i}$
$\displaystyle\hskip
11.38092pt\textbf{if}f(\mathbf{p}_{\text{best},i})<f(\mathbf{g}_{\text{best}})$
$\displaystyle\hskip
22.76228pt\mathbf{g}_{\text{best}}\leftarrow\mathbf{p}_{\text{best},i}.$

The use of $\mathbf{g}_{\text{best}}$ in the description of PSO above
represents the swarm _communication topology_ —in this case all members of the
swarm communicate and share knowledge of the global value
$\mathbf{g}_{\text{best}}$ and so this is referred to as the global
communication topology (hal_spso; engel). The most effective communication
topology was heavily researched during the development of PSO. Several other
communication topologies were considered. A random local topology, where
communication neighborhoods are randomly instantiated for particle $i$ with a
small number, say 3, of expected connections (i.e. particle $i$ can only share
the information its gained about $f$ to the particles in its neighborhood) was
found to be effective in improving the robustness of the swarm to entrapment
in local optima for multimodal objectives and exact optimal design problems
(engel; walsh). In this approach if the swarm does not find a better fitness
on $f$ at time $t$ then the communication links are randomly re-instantiated
at time $t+1$, thus creating a new local communication neighborhood for
particle $i$ and, hence, new information sharing among the particles. In this
scenario the term $\mathbf{l}_{\text{best},i}$, which represents the best
position in $\mathcal{X}$ found _in the communication neighborhood of particle
$i$_ to time $t$, replaces $\mathbf{g}_{\text{best}}$. Thus, the local
communication topology encourages the swarm to explore $\mathcal{X}$ more
thoroughly at the expense of more iterations. However, we view this as an
acceptable tradeoff if we are seeking the global optima of a difficult-to-
optimize objective function.

###  3.1 Extending PSO to the Optimal Design Problem

PSO is formulated to optimize a function $f$ which takes a $K\times 1$ vector
input $\mathbf{x}$. In contrast, optimal design requires optimizing an
objective function which takes an $N\times K$ dimensional matrix input
$\mathbf{X}$. We considered using the standard formulation of PSO by
implementing $\mathbf{x}^{*}=\text{vec}(\mathbf{X})$ where vec is the
vectorization operator (i.e. it stacks the columns of a matrix) as this would
give us access to off-the-shelf implementations of PSO. Due to programming and
research considerations we implemented an alternate but equivalent approach by
extending the PSO update equations as

$\displaystyle\mathbf{V}_{i}(t+1)$ $\displaystyle=\omega\mathbf{V_{i}}(t)$
(inertia) (13) $\displaystyle+$ $\displaystyle
c_{1}\mathbf{U}\odot\left(\mathbf{P}_{\text{best},i}-\mathbf{X}_{i}(t)\right)$
(cognitive) (14) $\displaystyle+$ $\displaystyle
c_{2}\mathbf{U}\odot\left(\mathbf{L}_{\text{best}_{i}}-\mathbf{X}_{i}(t)\right)$
(social) (15) $\displaystyle\mathbf{X}_{i}(t+1)$
$\displaystyle=\mathbf{X}_{i}(t)+\mathbf{V}_{i}(t+1)$ (position update) (16)

where the $N\times K$ matrix $\mathbf{U}=\\{u_{nk}\\}$ with
$u_{nk}\overset{i.i.d.}{\sim}U(0,1)$ (i.e. all elements are randomly drawn
from a uniform distribution). It is evident that equations (13-16) reduce to
Eq.s (9-12) when $N=1$ and the resulting vectors are transposed. Thus, in our
formulation, a particle is now a candidate design matrix $\mathbf{X}_{i}$ and
has a corresponding velocity matrix $\mathbf{V}_{i}$ which governs how
particle matrix $i$ will step through the search space
$\mathbf{X}\in\mathcal{X}^{N}$ in search of the exact globally optimal design.
To fully implement PSO to search for optimal designs on the hypercube, there
are several other technical problems to address (e.g. how do we deal with
design points that leave $\mathcal{X}$ during the search?). For these we have
adapted standard PSO solutions as found in the literature (hal_spso). Our full
PSO algorithm for generating optimal designs on the hypercube is provided in
Algorithm 1. The user can specify lower and upper bounds on each of the $K$
factors, and in line 1 we have set defaults as
$\mathbf{l}_{b}=-\mathbf{1}_{K}$ and $\mathbf{u}_{b}=\mathbf{1}_{K}$
respectively without loss of general study because the initialization of PSO
takes into account the size of the search space (a type of standardization),
but these may be relaxed to irregular hyper-rectangular regions if required
for alternate searches. Lines 4-11 show the initialization of the particle
swarm (which is here a set of randomly drawn design matrices and initial
velocities). In line 4 the $j$th row of initial design matrix is drawn
corresponding to a uniform distribution with lower and upper bounds
$\mathbf{l}_{b}$ and $\mathbf{u}_{b}$ (hal_spso). The $j$th row of the
velocity matrix of particle $i$ is initialized similarly as shown in line 5.
Line 6 shows that the stepsize of elements of the design matrices is
controlled for each of the $k$th factors via parameter $\nu_{k}^{\text{max}}$,
and hal_spso recommends setting this parameter to half the length of the
search space in one dimension though they do indicate that the use of this
parameter is superfluous with the use of the values of $\omega,c_{1},c_{2}$ as
employed in this paper (because the initialization of PSO takes into account
the size of the search space which is a type of standardization). Line 7
indicates the instantiation of $\mathbf{N}_{i}$ which represents a vector
indicating which particles are in the communication neighborhood of particle
$i$. Further, this vector is created by subroutine genNeighbors according to
hal_spso which sets neighborhood links for particle $i$ randomly such that the
expected number of communication links is 3 but may range from 1 (i.e. the
particle only informs itself) to $S$ (i.e. the particle can inform all other
particles). With initial particle positions and velocities determined, line 10
initializes the solution, and line 11 shows initialization of the neighborhood
best locations. Line 13 starts the swarm search which iterates until a
stopping criteria is met. Lines 14-16 check the solution at time $t$ and if it
is unchanged from time $t-1$ then the communication neighborhoods,
$\mathbf{N}_{i}$, are randomly re-instantiated. Lines 18-29 describe velocity
updating, position updating, and swarm knowledge updating as previously
discussed in this section. Note the subroutine confine on line 22. It is
possible that a particle moves outside the search space during the search and
so we implement a ‘reflecting wall’ particle confinement as in hal_spso.
Function confine checks every element of $\mathbf{X}$. If any element has
moved out of the search space, that position is set to the closest point on
the boundary (in its respective dimension) of the $\mathcal{X}$ and its
velocity is halved and reversed in that direction (i.e. it ‘bounces off the
wall’), see walsh for technical details. Last, lines 31-24 show how the local
neighborhood best locations are updated after the personal best’s
$\mathbf{P}_{\text{best},i}$ are updated. The algorithm progresses in this
manner until a stopping criterion is met.

1: Input: Objective function $f$, number of particles $S$, search space bounds
$\mathbf{l}_{b}=-\mathbf{1}_{K}$ and $\mathbf{u}_{b}=\mathbf{1}_{K}$

2: // Randomly draw $S$ candidate designs $\mathbf{X}_{i}$ and initialize

3: for each $i=1,\ldots,S$ do

4: $\mathbf{x}^{\prime}_{ij}\sim U_{K}(\mathbf{l}_{b},\mathbf{u}_{b})\text{
for }j=1,\ldots,N$ giving candidate design $\mathbf{X}_{i}$ with rows
$\mathbf{x}^{\prime}_{ij}$

5: $\mathbf{v}^{\prime}_{ij}\sim
U_{K}\left(\frac{\mathbf{l}_{b}-\mathbf{x}_{ij}}{2},\frac{\mathbf{u}_{b}-\mathbf{x}_{ij}}{2}\right)\text{
for }j=1,\ldots,N$ giving velocity matrix $\mathbf{V}_{i}$ with rows
$\mathbf{v}^{\prime}_{ij}$

6: $\\{v_{ijk}\leftarrow\min\\{v_{ijk},v_{k}^{\text{max}}\\}\\}$ for
$j=1,\ldots,N$ // limit stepsize

7: $\mathbf{N}_{i}\leftarrow\texttt{genNeighbors}(\mathbf{X}_{i})$ // generate
communication neighborhood for particle $i$

8: $\mathbf{P}_{\text{best},i}\leftarrow\mathbf{X}_{i}$ //set initial personal
best position

9: endfor

10:
$\mathbf{G}_{\text{best}}\leftarrow\underset{\mathbf{X}_{i}\in\\{\mathbf{X}_{1},\mathbf{X}_{2},\ldots,\mathbf{X}_{S}\\}}{\mathrm{argmin}}f(\mathbf{X}_{i})$
// current best design among swarm

11:
$\mathbf{L}_{\text{best},i}\leftarrow\underset{\mathbf{X}_{i}\in\mathbf{N}_{i}}{\mathrm{argmin}}f(\mathbf{X}_{i})$
for $i=1,\ldots,S$ // current best design in local neighborhoods

12: // swarm search loop

13: while stopping criteria not met do (this is iteration over $t$)

14: if $\mathbf{G}_{\text{best}}(t)=\mathbf{G}_{\text{best}}(t-1)$

15: // if the solution doesn’t improve, reform the local communication
networks

16: $\mathbf{N}_{i}\leftarrow\texttt{genNeighbors}(\mathbf{X}_{i})$

17: endif

18: for each $i=1,\ldots,S$ do // Update velocities and positions

19:
$\mathbf{V}_{i}\leftarrow\omega\mathbf{V_{i}}+c_{1}\mathbf{U}\odot\left(\mathbf{P}_{\text{best},i}-\mathbf{X}_{i}\right)+c_{2}\mathbf{U}\odot\left(\mathbf{L}_{\text{best},i}-\mathbf{X}_{i}\right)$

20: $\\{v_{ijk}\leftarrow\min\\{v_{ijk},v_{k}^{\text{max}}\\}\\}$ for
$j=1,\ldots,N$

21: $\mathbf{X}_{i}\leftarrow\mathbf{X}_{i}+\mathbf{V}_{i}$

22: $\mathbf{X}_{i}\leftarrow\texttt{confine}(\mathbf{X}_{i})$ // keep
candidates in searchspace

23: if $f(\mathbf{X}_{i})<f(\mathbf{P}_{\text{best},i})$ // update knowledge
about best known design to time $t$

24: $\mathbf{P}_{\text{best},i}\leftarrow\mathbf{X}_{i}$

25: if $f(\mathbf{P}_{\text{best},i})<f(\mathbf{G}_{\text{best}})$

26: $\mathbf{G}_{\text{best}}\leftarrow\mathbf{P}_{\text{best},i}$

27: endif

28: endif

29: endfor

30: for each $i=1,\ldots,S$ do // with personal best’s update, can update
local best

31: if $f(\mathbf{P}_{\text{best},i})<f(\mathbf{L}_{\text{best},i})$

32: $\mathbf{L}_{\text{best},i}\leftarrow\mathbf{P}_{\text{best},i}$ // update
best known design in local neighborhoods

33: endif

34: endfor

35: endwhile

36: Output: Particle swarm solution—the best optimal design found
$\mathbf{G}_{\text{best}}$

Algorithm 1 PSO for Generating Exact Optimal Designs on the Hypercube

###  3.2 Computing Architecture