A graph-theoretic framework for algorithmic design of experiments
Ben M. Parker Steven G Gilmour Vasiliki Koutra




 
1]Ben M. ParkerB.M.Parker@soton.ac.uk
2]Steven G Gilmour
1]Vasiliki Koutra
[1]Southampton Statistical Sciences Research Institute, University of 
Southampton, UK
[2]King's College, London

A graph-theoretic framework for algorithmic design of experiments
    [
    December 30, 2023
=================================================================


 In this paper, we demonstrate that considering experiments in a graph-theoretic manner allows us to exploit automorphisms of the graph to reduce the number of evaluations of candidate designs for those experiments, and thus find optimal designs faster.
 
 
We show that the use of automorphisms for reducing the number of evaluations required of an optimality criterion function is effective on designs where experimental units have a network structure. Moreover, we show that we can take block designs with no apparent network structure, such as one-way blocked experiments, row-column experiments, and crossover designs, and add block nodes to induce a network structure. Considering automorphisms can thus reduce the amount of time it takes to find optimal designs for a wide class of experiments.  

Keywords: Linear Network Effects Model, Optimal Design of Experiments, Automorphisms, Block Designs, Isomophic Designs
 


§ INTRODUCTION


In previous work <cit.>, a method for designing 
experiments on social networks was introduced. That paper presented a linear network effects model, which provided a framework for finding the optimal design of experiments when experimental units were connected according to some relationship, which was specified by an adjacency matrix, and treatment effects propagated to experimental units connected to the experimental unit to which the treatments were applied. 

Although in that paper the primary application was the design for experiments on social networks, a variety of examples were given showing how these networks could be useful in many applications, for example in agricultural experiments where experimental units were connected by some spatial relationship, and also in crossover trials, where experimental units were connected by temporal networks. 

In that paper, some methods were introduced to exploit symmetries in the network in order to improve the speed of finding an optimal design. In this current paper, we expand this argument, and show that there is a wide class of experiments that can be reformulated into a problem of design on the network, and that by presenting the problem in this networked form, we can suggest how to modify existing algorithms to allow optimal designs on the original problems to be found more quickly by exploiting automorphisms of the network. We recap the previous literature, and the paper we base this work on (<cit.>) in what remains of this section. In Section <ref> we give a brief overview of algorithms used currently for design, and propose a new algorithm for social networks based on the properties of the graph. We give examples of use of the new algorithm in Section <ref>, before showing how this algorithm is useful in a wide variety of design problems, even when there is no obvious network structure. We conclude briefly in Section <ref>.




 §.§ Previous literature


Due to the increasing prevalence of (particularly) social networks, the general field of network science, and statistical analysis of networks has shown growth over recent years. <cit.> provides a review of the importance of statistical inference and design in social networks, describing that although  the literature on analysis of networks has increased rapidly in recent years, there is limited research on experimental design on networks. Much of the research tends to relate to large-scale properties of networks, for example <cit.> draw inspiration from model-assisted survey sampling to provide a detailed view of how treatments might be chosen when there are structured relationships represented as networks between experimental units, specifically under a normal-sum model where the response of a node is governed by the sum of normally distributed random variables, determined by the responses of the node's neighbours. 


There have been some attempts to make use of graphs as a tool for finding optimal designs, mostly in the context of block designs. <cit.>, for the application of designing efficient experiments for dual-channel microarrays, presents a way of representing experiments with blocks of size two, where when a direct comparison is made between two treatments in a block it is represented as a link in a network. Designs (for the A_s-optimality criterion as we refer to it in this paper) were found by simulated annealing for up to n=20 nodes (treatments in the representation of <cit.>). <cit.> extends this work on microarrays, and presents optimal designs. <cit.> generalise this concept to more general block designs by considering the concurrence graph of a block design where nodes representing treatments are joined by a number of edges equal to the number of blocks in which the treatments appear together. By considering properties of the concurrence graph thus formed, they present some interesting combinatoric parallels, and methods for finding optimal designs for these blocked experiments. <cit.> present some parallels between optimal design criteria and well-known properties of graphs using concurrence graphs, in the context of block designs.

The idea of considering isomorphisms of designs - where two or more designs are equivalent- has been used to reduce the search space for finding optimal designs is established: for example <cit.> provide an algorithm for fractional factorial designs which makes use of isomorphisms to reduce the complexity of calculations; <cit.> provide a way of classifying orthogonal arrays by isomorphisms. <cit.> show that the problem of finding whether block designs are isomorphic is equivalent in complexity to the graph automorphism problem. Perhaps surprisingly, this result does not seem to have been used in a constructive method to find designs, as we do in this paper.


Most of the results in this present paper extend <cit.>, where we present a linear network effects model which describes the response of each node as depending on the treatment given directly to that node, and also depends on the treatments given to the neighbours of that node. This is summarised in Section <ref> below. Further recent work in this area follows in <cit.>.




 §.§ Recap of design for networks and notation

<cit.> considers networks of experimental units; a network G=(N,E) is an undirected graph, a collection of nodes N and edges E ⊆( N × N) where the nodes represent experimental units on each of which we apply some treatment. The edges represent relationships between the experimental units. 

We assume that if a relationship exists between two experimental units, the
response of each experimental unit is dependent on the treatment applied to the
other according to a linear network effects model which we specify below.

We have |N|=n experimental units, and m treatments. The
relationship between experimental units is
specified by the n× n adjacency matrix A where A_ik=1 if and only 
if i and k are
related and A_ik=0 otherwise. By convention, A_ii=0. 

We assume there is a “subject effect” on the response from experimental unit i of
τ_j
when
treatment j is given to that subject, and a “network effect” of
γ_l
if
treatment l is given to a connected experimental unit k (if A_ik=1). We
assume
that each experimental unit receives
exactly one treatment. 
 

We measure the response Y_i for each of our experimental units. Let t(i) be the
treatment applied to experimental unit i; then our response is modelled as

    Y_i= μ+τ_t(i) +∑_k=1^nA_ikγ_t(k)+ϵ_i.

We call this model the “Linear Network Effects Model” (LNEM). We assume that errors ϵ_i  are independent and
identically distributed with mean 0 and constant variance σ^2.
We assume that we wish to estimate the subject and/or network effects, or some
contrast of them.

We let ǔ_j be the indicator vector
with i-th element equal to 1 when treatment j is applied to subject i, and 
otherwise 0. In matrix form,

    (Y̌)   =[   1̌  ǔ_1  ǔ_2    …  ǔ_m Aǔ_1 Aǔ_2    … Aǔ_m ][  μ; τ̌; γ̌ ]=Fβ̌,

where F is our design matrix, and our
vector of parameters is

    β̌=(μ  τ̌^̌Ť  γ̌^̌Ť)^T=(μ  τ_  1 …τ_m  γ_1…γ_m)^T.







We calculate the Fisher information matrix as I=F^TF, and choose an optimality criterion: we seek to minimise the average variance 
of all pairwise differences of treatment
effects,

    2/m(m-1)∑_j=1^m-1∑_l=j+1
    ^m(τ_j-τ_l).
 This is defined to be A_s-optimality 
for
estimating the
differences in the treatment effects. To ensure estimability, without loss of generality we set τ_m=0.







§ ALGORITHM


We recap briefly some general properties of design algorithms, which allow us to put our algorithm in context.
Let the design space, the set of all possible assignments of treatments to experimental units for our experiment, be 𝒳. To assess how good a design x∈𝒳 is we calculate the value of some optimality criterion f(x), and seek to find the design(s) which maximise f(x); if our optimality criterion is such that we wish to minimise f(x) without loss of generality we maximise -f(x) instead. That is, we wish to find x^*=max f(x), which we call an optimal design. Typically f() will be some function of the Fisher information matrix, I(x) ; for example A-optimality, a popular criterion, involves taking the average variance of all unknown parameters, which corresponds to f(x)=(I^-1(x)), the trace of the inverse of the information matrix. D-optimality, which minimises the joint confidence region for the unknown parameters, can be expressed as f(x)=(I^-1(x)). Computationally, f(x) is normally hard to evaluate. For example, in A-optimality both the calculation of this information matrix and the resulting inversion of it will be computationally hard.

Typically when finding optimal designs, the design space 𝒳 which we must search over to find a design may be large. For example, if we have n experimental units which each take values in some set ℱ then the size of the design space will be |ℱ|^n, and we suffer from the “curse of dimensionality”. Even if F is the binary set {0,1} such that each experimental unit may be assigned one of two treatments, then the size of the design space may be 2^n. If F is a continuous set (e.g. ℛ or, a set of equivalent size, the interval [0,1], the design space is (infinitely) bigger.

Ideally, we would evaluate the optimality criterion f(x) at all possible designs x∈𝒳 to find the optimal design x^*=max_x∈𝒳f(x), but for large design spaces and/or criteria that are difficult to calculate, computational restrictions may mean we can not do this or at least not compute this in a reasonable time. In some cases, analytic results will enable designs to be found exactly without a search algorithm. However, in general we seek some algorithm that enables us to find the optimal design (or a design which is near optimal). In overview, existing algorithms for design seek to do one or both of the following: 


 
  * reduce the complexity of the calculation of the optimality criterion f(.)
 
  * evaluate a subset of 𝒳 such that the overall number of calculations of f(.) is smaller.


Many algorithms exist within these categories. The first category includes emulation, used extensively within computer designs, where instead of evaluating a complicated function f(), we evaluate an emulator g() , a function which is simpler to evaluate but we believe has similar properties such that f(x)≈ g(x) for all x ∈𝒳. Updating formula can also be used in some problems, where we know that f(x+z)=f(x)+h(z) for some easy-to-calculate function h. 

The second category of algorithms include space-filling designs, where a representative sample of designs that “fill” the space  are evaluated. There are also many optimisation algorithms, which use the t previously evaluated designs x_0,x_1,… x_t-1 and their optimality criterion function evaluations f(x_0),f(x_1),…,f(x_t-1), to choose which design x_t should be evaluated next.  Fedorov exchange algorithms, coordinate exchange algorithms, simulated annealing, particle swarm optimisation, many stochastic search algorithms such as Nelder-Mead, are the subject of research. Some of these algorithms are deterministic, in that the choice of x_t is mandated, others are stochastic in that x_t is chosen randomly. To guard against designs which are local maxima, the initial x_0 is often chosen randomly, even if the rest of the algorithm is deterministic, and the algorithm repeated with many random x_0 chosen as starting designs.

Much of the literature in experimental design revolves around trying to find better algorithms for particular problems. In this work, we do not seek to find new algorithms, but to know how mapping the optimisation problem to a network domain allows us to vastly reduce the number of designs evaluated in order to allow us to find better designs. We can use many existing algorithms with small modifications in this new network framework.



 §.§ Using networks for reducing the state space

We show in this section how we can use the network representation of the problem as an advantage when we consider algorithms for design.  <cit.> introduced some ideas of how we might improve our algorithms, and we recap these here briefly before expanding these ideas and presenting a new algorithm. 



  §.§.§ Symmetry of labels

For many optimality criteria we are only interested in differences between treatments, the 
treatment effects themselves are irrelevant, and treatments are equivalent up to
relabelling. For example, designs {1,2,2,3,1} and {1,3,3,2,1} are equivalent. We can
thus reorder any design so that we only evaluate designs where the first
occurrence of label j must come before the first occurrence of label j+1.
Without loss of generality, we can assign treatment 1 to experimental unit 1.
 
 



  §.§.§ Symmetry of networks





In Figure <ref>, subgraphs 1 and 2 are exchangeable; i.e. if
we consider subdesign A on subgraph 1, and subdesign B on subgraph 2 (call this
[A1,B2]), we need not also consider [A2,B1] as by symmetry this design has the
same criterion function value. 

We can thus reduce our design space greatly if we can identify subgraphs where
the designs are exchangeable. This is equivalent to finding an
automorphism of our network, a relabelling or permutation of the set N
such that the edges E are preserved.  This is known as the Graph Automorphism 
Problem, and finding efficient algorithms to solve this problem has been the topic of much research <cit.>. Until recently, the general case of this problem was thought to be not solvable in polynomial time, that is no algorithm existed that could find a solution for a graph with n nodes in O(n^k) time for some constant k. (Recently, an algorithm has been found that claims a solution exists in quasi-polynomial  (exp((log n)^O(1))) time (<cit.>, although this result has not yet been peer-reviewed.). 

Whatever the theoretical complexity, in practice fast algorithms exist already that can effectively and quickly find automorphisms in most cases, for example the VF2 algorithm<cit.>, which in the general case finds isomorphisms between two graphs G_1 and G_2. This algorithm is based on a tree search, where a set of nodes corresponding to a partial match between subgraphs of G_1 and G_2 is maintained, and then nodes connected to these subgraphs are considered to see if they can extend the matching subgraphs. This algorithm is fast, and implemented by the popular  package (<http://igraph.org/>) available in many programming languages, including R. Here, we find automorphisms by setting G_1=G_2. Many algorithms exist, but we chose VF2 as it was fast enough for our purposes, and was readily incorporated in the software made available with tis paper.



 §.§ New framework for design algorithm


We present below a general algorithm:


  * Rewrite original problem in graph form.



  * Find the automorphisms  of the graph using the VF2 algorithm. This is a set such that for any labelling of the graph corresponding to design x, all designs isos(x) will be automorphisms of x. 

  * Set k=1, num_=0 and pick an initial candidate design x_1. 


  *  Check whether x_k=min_(x_k), i.e. if x_k is the first design in the lexicographical order of automorphisms of x_k. If so, 
 
 
  * evaluate the optimality criterion d_k=f(x_k)
 
  * increment num_ by 1. 



  * Set k=k+1. Pick next candidate design x_k according to algorithm of choice  based on x_1,…,x_k-1 and d_1,…,d_k-1, so that x_k=[(x_1,…,x_k-1),(d_1,…,d_k-1)]. 

  * If we have reached a stopping criterion, i.e if stop)[(x_1,…,x_k),(d_1,…,d_k),num_]=1 for some stopping criterion , then stop, otherwise repeat from <ref>.



Typical choices of the  algorithm might be the Fedorov exchange algorithm, an exhaustive search, etc.
Typical choices of the  criterion might be that k=k_1 or num_=k_2 representing some fixed budget on the number of design points examined or function evaluations made, or that max (d_k-t,…,d_k)=d_k-t, i.e. we have not seen an improvement in the last t steps.


  §.§.§ Treatment Structure

In this paper we restrict ourselves to designs on m unstructured treatments; that is our design space is ={(x_1,x_2,…, x_n| x_i∈{1,…,m}∀ i} such that ||=n^m. We do this for clarity of exposition of our method, but the general ideas extend to factorial and other treatment structures.




  §.§.§ Illustrative example

Consider the network represented by 3-1-2, where we have three nodes such that nodes 3 and 2 are each connected to node 1 and there are no other connections. It is clear that the network is automorphic to 2-1-3. Let 𝒯={A,B}, and let a design be represented by x_1x_2x_3 where x_i∈𝒯 is the treatment given to node i. In other words ABB means we give treatment A to node 1, B to node 2, and B to node 3. Note that 𝒳={AAA,AAB,ABA,ABB,BAA,BAB,BBA,BBB} where the elements of the design space are clearly in lexicographical order. Let, X_1=AAA, and we let  be the next design lexicographically, and let our stopping criterion be  if x_k=BBB.  Then the automorphisms allow us to show that two pairs of designs are equivalent (AAB=ABA, BAB=BBA); thus in step <ref>, we do not evaluate f(x) for designs ABA or BBA. We have therefore reduced the number of designs we must consider from 8 to 6..



  §.§.§ Rationale

For many designs, evaluating the optimality criterion function f() requires (sometimes substantial) computational time, and by reducing the number of evaluations required, we can reduce the computational time required. Alternatively, for the same computational time, we can search more of the design space and hope to find a better design. 

There is some computational overhead in the new algorithm; evaluating the automorphisms initially can be slow, however this needs to be done just once for each network, so for large designs this time is not a large proportion of the total computation time. There is also some overhead in step 3(a), as checking whether a design is lowest in lexicographical ordering of all automorphisms of that design requires some computation. However, we show that this overhead is substantially lower than the cost of evaluating f() in some simple examples, so our algorithm can lead to substantial improvements in finding designs. For particular networks, for example networks with a large number of automorphisms relative to the size of the network, the overhead may be prohibitive. This may be the case for dense or highly regular networks. We investigate the computational complexity of the algorithm in <ref>.





§ EXAMPLES OF USE OF AUTOMORPHISMS IN FINDING DESIGNS



We return to a small example which was Example 1 in <cit.> and represents a social network. This is shown in Figure <ref>, together with the optimal design for estimating with minimal variance the difference between the m=2 treatment effects, τ_1-τ_2, under the LNEM (<ref>). The list of edges is shown in the Appendix.

The optimal design was evaluated using exhaustive search. For each candidate design in our design space (there are ||  =2^10 designs), the information matrix must be calculated (this a 4 by 4 matrix with columns representing μ, τ_1, γ_1, γ_2, recalling we can omit the column corresponding
to τ_2=0). This matrix must be inverted, and the inversion is responsible for much of the computational effort in finding the optimal design. In practice, we can note that as we are only interested in differences between subject effects, we can appeal to symmetry of labels to fix the treatment for experimental unit 1 as the first treatment.

In <cit.> we did not use the symmetries of the network, but found the design under exhaustive search; here we use an exhaustive search (so that our  design calculates the lexicographically subsequent design), but will only evaluate designs which are first lexicographically amongst all designs which are automorphisms of the current considered design. We will thus evaluate fewer designs, and  when one design from each automorphism class has been evaluated.

We naturally find the same optimal design, but evaluate the information matrix 236 times as opposed to 507 times, and use a processing time of 0.02 seconds as opposed to 0.04 seconds.




We can repeat this evaluation for all the networks in <cit.>, the edge list for each of which is shown in the Appendix. We chose these networks as they represent a variety of different types of network structures, which relate to a variety of different applications. The results are shown in Table <ref>, and indicate that the number of function evaluations is greatly reduced for all the networks (except network 2, where there are no non-trivial automorphisms) and that the time taken to perform the exhaustive search is greatly reduced, even when taking into account the time needed for calculating automorphisms and lexicographical ordering checking as described above. 

We include the computational time as a guide; this will vary greatly between implementations in different software for these small networks, but gives an idea of how effective computationally the use of automorphisms is for finding designs on network problems, even for small networks.




 §.§ Implementing a simple coordinate exchange algorithm


We demonstrate how the framework can be combined with a non-trivial design algorithm by showing as an example the cyclic coordinate descent algorithm described fully in<cit.>.

We pick a random starting design. Each node is considered in turn, and for each node we cycle through all possible treatments for that node to see if changing the treatment allocated to that node can improve the design. If we find an improvement, we fix that design, and start again. If after checking all treatments we find that changing the treatment on this node does not improve the design, we move to the next node. When we have cycled through all nodes without improvement, we stop.

In the language of our general algorithm above, we pick

    x_k   =   [(x_1,…,x_k-1),(d_1,…,d_k-1)]
       =   [x^*_i+u_k-i,m] m
 
where u_k,m is the unit vector with all elements zero except the (k/m)-th, which is 1, and x^*_i is the maximum design so far discovered such that x_i^*=max_1≤ i ≤ k-1f(x_i). for design i. The stopping criterion is chosen such that  


    [(x_1,…,x_k),(d_1,…,d_k),num_]=1  k>nm  d_k=d_k-nm,
 i.e. if we have evaluated all designs resulting from changing one coordinate and not seen an improvement.

 
 
 

We perform 100 random starts to mitigate the algorithm finding local maxima. To demonstrate how this algorithm can work with automorphisms, we apply this algorithm to the same 6 networks as in Table <ref> above. We present the results as Table <ref>. We can see that even this simple algorithm is compatible with the network framework of taking automophisms into account, and a very small number of function evaluations are made, with generally highly efficient designs being found.


 

§ BLOCK DESIGNS AS NETWORKS


We seek to generalise the framework developed to include a wider class of models. We start with a simple example.

Suppose we have four blocks, each of four experimental units. Let us suppose we wish to allocate treatments labelled {1,…,m} with corresponding unknown effects {τ_1,…,τ_m} to the experimental units to maximise some optimality criterion; for example we may wish to estimate 

    2/m(m-1)∑_j=1^m-1∑_l=j+1
    ^m(τ_j-τ_l).
 This is defined to be A_s-optimality for
estimating the
differences in the treatment effects. 
We call this criterion function ϕ_1.

A typical schematic for the structure of the experiment might be written down as in the left of Figure <ref>. It is clear that experimental units 1,2,3, and 4 are contained within block 1, and so on. We might choose to represent the experiment pictorially as in the right of Figure <ref>.





Why is this useful? Let us now imagine that the 20 nodes {1,…,16,B1,B2,B3,B4} are experimental units that are connected according to the relationship shown in Figure <ref>. We forget that the nodes have any special meaning now, save that we may assign treatments 1,…,m only to the first 16 nodes, and that node B1 is assigned treatment m+1, node B2 assigned treatment m+2, B3 assigned m+3 and finally B4 assigned m+4.

In order to find the A-optimal design, we can write the model for our blocked experiment as  

    Y_i= μ+τ_t(i) +∑_k={1,…,16,B1,B2,B3,B4}A_ikγ_t(k)+ϵ_i,    i=1,…,16

where A is the adjacency matrix where A_ij=1 whenever there is a link as shown in the Figure (e.g. A_14,B4=1, and A_ij=0 otherwise).

By writing b_j(i) as ∑_k={1,…,16,B1,B2,B3,B4}A_ikγ_t(k), we can see immediately that this is equivalent to

    Y_i= μ+τ_t(i) +b_j(i)+ϵ_i,
 a more familiar representation of a blocked experiment, where b_j(i) is block effect of experimental unit i in block j. As A_ij=1 if and only if node i is linked to node j, and this only happens when experimental unit i is in block j for j={B1,B2,B3,B4}, we can  replace γ_t(B1)=γ_m+1=b_1, γ_t(B2)=γ_m+2=b_2, and so on. 

Perhaps more simply, the network effects of the special nodes representing blocks in the linear network effects model are replaced by the block effects.  We can think of the block effects as propagating like network effects from the special nodes B1, B2, B3, B4. We will normally wish to estimate some function of the treatment effects only (block effects are often not of interest but it is important to account for them in the design), and as the treatment effects for nodes B1, B2, B3, B4 are irrelevant (we cannot measure the experiment units B1, B2, B3, B4 directly as they are not real units!), we ignore them and use the same optimality criterion as before. As we are not estimating all effects, strictly this is now A_S optimality rather than A optimality. 


The properties of the two designs can be summarised in Table <ref> .



We have thus shown (at least for this simple example) an equivalence between a traditional design (in this case a block design) and a network design as found in the paper <cit.>. A major advantage of writing a block design in this way is that we can now find optimal designs for block designs in the same way as for network designs, and that we can use the automorphisms found from the network representation of the experimental framework in order to find optimal designs faster. The network setting allows us to develop one set of effective algorithms, rather than regarding blocked experiments as a special class  with its own method for finding optimal designs.



 §.§ Other block designs as networks


 

In a similar method to the one-way block design, it is possible to represent a variety of block designs as networks. For all blocking factors each block can be written as a new network node, linked to all experimental units within that block. We present examples for a double-blocked experiment (a row-column design) as Figure <ref>. This is readily extendable to more than two blocking factors.





In a crossover design, subjects receive treatments sequentially over several time periods, and an experimenter wishes to account for the treatment possibly still having some effect in a later time period than that for which it is applied. We can represent a crossover-design similarly to a row-column design as Figure <ref>. Experimental units are subject-period combinations, and we assume that subjects behave similarly, so subject is a blocking factor, and periods may be a blocking factor, so we add extra block nodes as we did for the row-column designs. In addition, we assume a treatment applied to a subject may affect the subject in the next time period, so there may be a carryover effect from experimental unit 1 to 2, 2 to 3, 4 to 5, etc. We represent this carryover effect by the network effect of the treatment given to the experimental unit preceding in time, but note this carryover effect only extends to the next time period, so links in the network are directed. This is similar to Example 6 in <cit.>.





 §.§ Advantages of changing design space

The argument explained in <ref> applies: by evaluating only one design from each class of automorphic designs, we can make significant savings in the amount of time to evaluate candidate design. 

We evaluate several experimental structures to demonstrate this in practice, for which the experimental designs are well-known:

 
  * 3 blocks of size 3 (n=9), with 3 treatments. (The optimal designs are randomised complete block designs.)
 
  * 4 blocks of size 3 (n=12), with i) 3 and ii) 4 treatments. (The optimal designs are i) randomised complete block designs and ii) balanced incomplete block designs.)
 
  * A row-column structure with 3 rows and 3 columns, each row-column intersection containing a single experimental unit, with 3 treatments. (The optimal designs are Latin Squares of size 3.)
 
  * A row-column structure with 4 rows and 4 columns, each row-column intersection containing a single experimental unit, with i)3 and ii)4 treatments. (The optimal designs in ii) are Latin Squares of size 4.)




We do not claim that these designs would be sensibly found via this method, as the solutions are known analytically, but we seek to demonstrate the benefits of using automorphisms in reducing computational time.

The results are presented as Table <ref>. Clearly the number of optimality function evaluations is vastly reduced for networks a large number of automorphisms. However, for some networks, such as 2i and 2ii, the overhead in the implementation caused by searching for automorphisms and then checking the lexicographical order of the candidate designs might make the use of automorphisms inefficient. 

For larger designs, such as 4i, with a moderate number of automorphisms, we see a more than tenfold reduction in processing time and a hundred-fold reduction in number of evaluations, suggesting that this method works better for networks with a moderate number of automorphisms. 

In practice, by taking into account network structure, practitioners may be able to make a very large saving in time to find optimal designs with very little modification to existing algorithms, meaning designs found using stochastic algorithms in fixed computing time may be better. 





 §.§ Computational Justification


Recall that ||=n^m for unstructured treatments. For our A_s optimality criterion, and indeed most common optimality criteria, when calculating f(x) we must calculate the Fisher information matrix and invert it to find a variance-covariance matrix, then make some calculation of this final matrix. We assess the typical computational complexity for each design we consider.

 
  * The complexity for calculation of the Fisher information matrix from F^TF where X is an n×(2m) matrix is typically O(n^2m).
 
  * The Fisher information matrix in our model (<ref>) is of size 2m×2m. Inverting the matrix depends on the algorithm used, but is O(m^k) where 2<k≤3.
 
  * Calculating the optimality criterion from the 2m×2m variance-covariance matrix depends on what criterion is used; the trace (A-optimality) is O(m), taking the determinant (D-optimality) will typically be O(m^k) where 2<k≤3.
 

  Thus for each design evaluated, and assuming a simple optimality criterion, as n>>m (in general we have many more experimental units than treatments) the limiting step is the first step above and therefore calculating f(x) has computational complexity of O(n^2m).
  
  The complexity of the overhead in the new framework algorithm involves i) the initial time to calculate automorphisms originally, and ii) the computational cost of checking whether each design is lowest lexicographically amongst all possible automorphic designs. 
  
  Let us assume the size of  is z. i.e. there are z automorphisms for our network. We map each design of length n to an automorphism of itself which reorders the design, each time this operation is O(n). We must do this z times (once for each automorphism), and then sort the resulting list of automorphic designs to find the smallest, which can be done in O(zlog z) for a good sorting algorithm. Thus the overall computational complexity of checking whether a design is first lexicographically is O(nz + zlog z).


  Thus the computational complexity involved in the proposed algorithm where we first check that a design is valid, and if it is calculate the design, is O(nz + zlog z+ n^2m/z), where the z in the denominator in the second term arises because we must calculate f(x) for one design in every z.
  
  Evaluating one design for the algorithm without automorphisms is O(n^2m); for the new method is O(nz+zlog z+ n^2m/z]). The ratio of the new to the old is 
    O(z/nm+zlog z/nm+1/z)≈ O(zlog z/nm+1/z)

  
  Thus the effectiveness of our new algorithm depends on the relative sizes of n, z, and m. We find z must be small compared to nm but not too small. This is supported by the results in Table <ref>, which show the computational time to be reduced significantly for moderate z, but actually increased for very large z, as the overhead in finding and ordering automorphisms is high here.
  













§ CONCLUSIONS


We have shown that the use of automorphisms for reducing the number of evaluations required of an optimality criterion function is effective on designs where experimental units have a network structure. Moreover, we have shown that we can take block designs with no apparent network structure, such as one-way blocked experiments, row-column experiments, and crossover designs, and add block nodes to induce a network structure. Using automorphisms on these experiments with induced networks is also effective at reducing the complexity of experimental design algorithms.


From a practical point of view, many algorithms for design exist in isolation; we must program an algorithm for split-plot designs differently, perhaps, to how we program row-column designs. We argue that a framework such as we suggest with this paper may promote general purpose algorithms, neggating the need to maintain different algorithms for particular designs. Although algorithms for design can take into account automorphic designs to avoid recalculation of f(x) for an automorphic design, in general they do not, and this network representation allows automorphisms to be found readily by algorithms that are quick and available in existing software. 
We believe that this work may lead to standardisation of algorithms across seemingly different classes of experiment.

Harvard




§ APPENDIX


The following lists the edges for the six examples shown in this paper. These are also shown visually in examples in <cit.>. “” means an edge exists between i and j, and the corresponding entry in the adjacency matrix A(i,j)=A(j,i)=1. If no edge it shown then A(i,j)=A(j,i)=0. Examples 1 to 5 are undirected networks. In Examples 6, the network is directed, such that “” means that A(i,j)=1, but note that A(j,i)1.





§ SOFTWARE


A draft package to find designs for networks, including the augmented networks for block designs, is available at <https://www.dropbox.com/s/fgfn2f2my4yqlzr/networkDesign_0.0.0.9001.tar.gz?dl=0>. Scripts are provided to allow the reader to reproduce many of the results in this paper, as well as those in <cit.>.


