Generating designs for comparative experiments with two blocking factors

Often, comparative experiments involve a single treatment factor and two blocking factors, for example, augmented row–column, two‐phase, and incomplete row–column experiments. These experiments are widely used in agriculture. Finding good designs for these experiments is a major challenge when the number of treatments is large and the blocking structure is complex. In this paper, we first propose a new search algorithm that is combined with efficient update formulae, so that optimal designs with two blocking factors can be found within a reasonable time. Second, we compare augmented row–column designs generated with our new method to those obtained from CycDesigN, DiGGer, and the OPTEX procedure of SAS in terms of computing times as well as the quality of solutions. Third, we illustrate our proposed approach with four applications. We show an example where our efficient update formulae work while existing update formulae cannot be applied, and we use our search framework to generate augmented row–column, two‐phase, and incomplete row–column designs. We end the paper with a conclusion along with suggestions for potential applications.


INTRODUCTION
Many comparative experiments in agriculture, pharmaceutical industry, and engineering are very expensive and time-consuming (John & Williams, 1995;Bailey, 2008).Thus, it is crucial to use a good design that allows collecting high-quality data with minimal correlation (Dean et al., 2015).To set the stage, consider the following three examples from plant breeding that motivate our work: Example I (augmented row-column design): At an early stage of a plant improvement program, new varieties (i.e., treatments) can only be tested on one experimental unit due to limited availability of seeds, while standard varieties can be replicated more often.The main goal is to test the new varieties.A design for such an experiment is called augmented design (Federer & Raghavarao, 1975;Federer et al., 1975).For example, a plant breeder may wish to study 21 varieties in a field experiment, where 18 new varieties, called entries, are unreplicated, and 3 standard varieties, called checks, are replicated 6 times each.The field layout has 6 rows and 6 columns.Since there might be row-to-row and column-to-column variation, it F I G U R E 1 An experimental layout for Example III with  1 = 11 rows,  2 = 16 columns, and 32 unusable plots (X).
is important to find a good arrangement of these varieties in the field layout taking the row and column effects into account.In experimental design jargon, the row and column factors are called blocking factors.
Example II (two-phase design): Often, agricultural experiments are run in two phases.Such experiments are referred to as two-phase experiments (McIntyre, 1955).Commonly, Phase 1 is performed in a field, while Phase 2 is performed in a lab.Consider an experiment in plant breeding.In Phase 1, 10 varieties are tested in a field trial involving 10 blocks of size 6.In Phase 2, 60 samples obtained from the experimental units in Phase 1 are analyzed in a lab.Each day we can only test four samples.Thus, 15 days are needed to test 60 samples.Since there might be field-block-to-field-block and day-to-day variation, it is crucial to take these sources of error into account as blocking factors.Unlike Example I, where the blocking structure is crossed, the ideal blocking structure of Example II is unknown.Thus, finding a good allocation of treatments to the blocking structures is challenging.It is therefore useful to study the merit of specific blocking structures for this example so that the complexity of this design problem can be reduced.
Example III (incomplete row-column design): Occasionally, in a field trial, for various reasons, it might happen that breeders perform an experiment in a rowcolumn layout where a few experimental plots cannot be used.We refer to a design for such an experiment as an incomplete row-column design.For example, Johnstone (2003, see Example 4 where the role of rows and columns is swapped) considered an experiment laid out in an 11row by 16-column array, where the preparation of some internal plots failed and rows of the layout had a gradation of dampness and shade which left some lower boundary plots unusable.There were 32 unusable plots (i.e., missing plots) in this row-column layout, which are indicated by "x" in Figure 1.The main goal of this experiment is to arrange 8 varieties, each of which has 18 replications, in this incomplete row-column structure so that we can estimate treatment effects precisely, correcting for the row and column effects.
Examples I, II, and III involve a single treatment factor and two blocking factors.In what follows, we use the term designs with two blocking factors to refer to designs for such experiments.Often, these experiments are performed several times in different trials and require substantial time to be completed.Therefore, it is crucial to use a good design allowing an insightful conclusion without many trials, so that time can be saved.Moreover, it is necessary to guarantee that all treatment effects can be estimated under a suitable linear model (see Section 2) and are confounded with the block effects as little as possible.Here, we focus on finding good designs with two blocking factors that possess these key features.
Generally, designs with two blocking factors can be viewed as generalized forms of row-column designs, where the rows (columns) correspond to levels of the first (second) blocking factor.Any cell intersected by a row and a column then might have zero, or more than one experimental unit.In the current literature, many studies focus on designs with two blocking factors with a crossed or nested blocking structure in which all cells are filled with one experimental unit and where there is equal replication for all treatments, see John and Williams (1995, Chap. 5).However, there are only few studies on computer search algorithms to tackle designs with two blocking factors and a complex blocking pattern (e.g., an incomplete row-column pattern).Particularly, it is unclear whether the existing computer search algorithms are indeed efficient and fast enough when searching for good designs in complex scenarios like the three sketched above, especially when the number of treatments is large, which is often the case in comparative experiments.Furthermore, given a known blocking structure, computer search algorithms typically perform many exchanges of two different treatments and, following each exchange, evaluate the Moore-Penrose inverse of the information matrix (Section 2).It is computationally expensive to find this inverse using standard procedures.Hence, it is desirable to derive efficient update formulae, and an efficient solution to this problem was given in John (2001).However, there are two cases in which this solution does not apply, as will be shown in this paper.Additionally, even though John's update formula has been used extensively (John et al., 2004;Nguyen & Williams, 2006;Williams & John, 2007), there are few studies quantifying the benefit of this update formula in terms of computing time.Moreover, it might be desirable to efficiently implement an update formula multiple times in parallel, for example, when many designs need to be evaluated simultaneously in searching for optimal designs in which blocking structures are given.All these issues are studied in this paper, whose main goal is to provide a unified search framework that can handle designs with two blocking factors, a complex blocking pattern, and equal or unequal replications within a reasonable computing time.
Section 2 proposes a linear model for designs with two blocking factors and criteria for selecting good designs.Section 3 proposes a new search framework for obtaining good designs with two blocking factors.Section 4 proposes an efficient update method that works also in cases where the method of John (2001) is not applicable.Section 5 evaluates the performance of our new framework for augmented row-column designs.Furthermore, in this section, we compare our designs with those generated by CycDe-sigN (VSNi, 2020), DiGGer (Coombes, 2020), and OPTEX (Pereira & Tobias, 2015;Piepho, 2015).Section 6 provides four further applications.Section 7 ends the paper with a conclusion and a suggestion for possible applications of our work.

LINEAR MODEL AND DESIGN CRITERIA
Suppose that we wish to find a good arrangement of  treatments in a known blocking structure involving two blocking factors and  observations.The th blocking factor has   levels ( = 1, 2).Treatment replications can be equal or unequal.A linear model for this design problem is of the form: where  is an  × 1 response vector,  is an  ×  binary design matrix,  is a  × 1 vector of fixed treatment effects,  = [ 1  2 ], where   is an  ×   binary design matrix for the th blocking factor,  = [ ⊤ 1  ⊤ 2 ] ⊤ , where   is the vector of fixed block effects, corresponding to the th blocking factor, and  is an  × 1 vector of random errors.We assume that the expected value of  is zero and its covariance matrix is  2   , where  2 is a constant variance and   is an identity matrix of size .Blocks will often be modeled as random at the analysis stage for the recovery of information.However, to consider blocks as random at the design stage, we would need to assign values to the variance components, unless it can be established that a design is optimal for all non-negative values of the variance components.As the variances are commonly unknown, it is one option to take a conservative approach and to design for the worst case, i.e., that which leads to the largest variances of estimated treatment effects.This occurs when the between-block variances tend to infinity, which is equivalent to the fixed blocks case (Bueno Filho & Gilmour, 2003).Thus, here, we use a linear model with fixed block effects to generate designs with two blocking factors.This model is also called a linear fixed-effects model.Also, without loss of generality, we assume  2 = 1 during the search of optimal designs.Then, the treatment information matrix () can be computed as where  =   − ( ⊤ ) −  ⊤ , and ( ⊤ ) − is a generalized inverse of  ⊤ .Here, we limit our attention to connected designs for which rank() =  − 1.

Design criteria
In comparative experiments, designs are often ranked according to an overall efficiency factor.For several classes of designs evaluated using linear models with fixed block effects, there is a link between the overall efficiency factor and the -optimality criterion, computed as the sum of the inverses of the non-zero eigenvalues of the treatment information matrix (John & Williams, 1995;Morgan & Stallings, 2014).If this link exists, designs obtained by maximizing this overall efficiency factor are -optimal.
We rank designs according to the overall efficiency factor

𝑟𝑉
, where  = ∕, and  is the average variance of all estimated pairwise comparisons.Designs with high  values are preferable. is given by where Ω is the Moore-Penrose inverse of , with entries denoted by   for 1 ≤ ,  ≤ ,   is a vector of ones of size , and  ⊤  Ω  = 0 since  ⊤    = 0 (John & Williams, 1995).Furthermore, it is straightforward to show that  is also equal to the harmonic mean of the non-zero eigenvalues of 1   (see Equation ( 3)).In what follows,  is also called the design efficiency.
For the purpose of comparisons with designs generated in CycDesigN, we use one of the  series of criteria (Williams & Piepho, 2015,   ), which we refer to as an overall weighted efficiency factor ( Ā).This criterion has weights proportional to treatment replications and is called the average efficiency factor in CycDesigN in the nonresolvable design menu with unequal replications.Designs with high Ā are preferable.The main difference between  and Ā is that the latter uses replication numbers as weights (Williams & Piepho, 2015, Equation (3)).As shown in Equation ( 2), the calculation of  requires a generalized inverse of  ⊤ , which is rather complex to compute.Hence, it is useful to derive a simple form of a generalized inverse of  ⊤  so that  is easy to compute.This will be useful when deriving our update formulae in Section 2. In the next two sections, we will derive simple forms for the information matrices of augmented row-column and two-phase designs considered in this paper.For incomplete row-column designs, it is infeasible to find a simple form of generalized inverse of  ⊤ .In this case, we use the Moore-Penrose inverse of  ⊤  to compute .

Information matrix of augmented row-column designs
For augmented row-column designs (ARCDs), we assume that there are   treatments with   = ( 1 ,  2 , … ,    ) replications, while  −   treatments have only one replication.The treatment information matrix for an ARCD is computed according to a modification of Equation ( 2).Here, we use where ⊗ denotes a Kronecker product.With these specifications, ]), where diagm() is a diagonal matrix with all main diagonal entries equal to .Note that it is possible to drop   in  * .However, the advantage of including it in  * is that ( * ⊤  * ) − can be expressed in the above convenient form.Then,  can be expressed as: where  = (   ⊤ −  ) ⊤ is the replication vector for the treatments, and  1 =  ⊤  1 and  2 =  ⊤  2 are row and column incidence matrices (Bailey, 2008).

Information matrix of two-phase designs
Unlike ARCDs, for two-phase designs, we assume that treatments are equally replicated ( replications).In Phase 1, we need to find a suitable arrangement of  treatments in  1 blocks of size  1 .In Phase 2, the  =  samples obtained from Phase 1 are then arranged in  2 blocks of size  2 .As shown in Vo-Thanh et al. (2023), two-phase designs can be constructed from nested rowcolumn designs.We describe a nested row-column design as follows.Let  be the greatest common divisor of  1 and  2 , and  be the greatest common divisor of  1 and From this, we have  =  1  1 =  2  2 .This implies that  1 =  2 and  1 =  2 .Based on these facts, we define  superblocks where each superblock is a row-column layout with  1 rows and  2 columns.In each superblock, the intersection of a row and a column is called a superplot and contains  experimental units ( ≥ 1).Thus, experimental units in a row form a block in Phase 1, while experimental units in a column form a block in Phase 2. We will illustrate the use of these parameters to construct a two-phase design in Section 6.2.The design matrices  1 and  2 then are (5) Remarkably, under this nested row-column structure, two-phase designs can be viewed as row-column designs.Unlike ARCDs, where every cell has one experimental unit, cells in a row-column design for two-phase experiments have either 0 or  experimental units.The following proposition provides a generalized inverse of  ⊤  for a nested row-column design: Proposition 1.Given  1 and  2 in Equation ( 5), we have where   is a zero vector of size  and  , is an all-ones matrix of size  × .
The proof of this proposition can be obtained easily by observing that  ⊤ ( ⊤ ) −  ⊤  =  ⊤ .When considering two-phase designs under the nested row-column structure, by combining Equations ( 2) and ( 6), the treatment information matrix can be expressed as: where   is an  ×  matrix of ones.
CycDesigN uses SA and focuses on linear fixed-effects models.The DiGGer (Coombes, 2020) and od (Butler, 2020) packages are based on RTS1 and RTS2, respectively, and focus on linear mixed models, where treatment effects or block effects are considered as random.The search algorithms used belong to the class of local search algorithms, some of which require many tuning parameters and have complex acceptance strategies.We refer to the textbook of Hoos and Stützle (2004) for those who are unfamiliar with SA and TS algorithms.Choosing appropriate values for the tuning parameters is often one of the major challenges when using these algorithms.Also, it is unclear whether these algorithms are efficient and fast enough when searching for ARCDs.Moreover, these packages are not suitable for finding optimal two-phase and incomplete row-column designs under the fixed-effects model considered in this paper.By contrast, the hill climbing (HC) algorithm, one of the earliest search techniques, does not require tuning parameters (Appleby et al., 1961).The main problem with the HC algorithm is that it does not accept inferior solutions.Also, it converges quickly to a local optimum, but the quality of its solutions is often not high (Hoos & Stützle, 2004).Thus, the HC algorithm aims at extraordinary intensification and suffers from inadequate diversification.Intensification refers to search algorithms that aim to greatly improve a solution, while diversification refers to search algorithms that help achieve a reasonable coverage in the search space (Hoos & Stützle, 2004 version of the well-known late-acceptance hill climbing search (Burke & Bykov, 2017, LAHC).Both algorithms aim to overcome the main downside of hill climbing.There are two advantages when using DLAS or LAHC.First, these algorithms require simpler tuning parameters, compared to the SA, TS, RTS1, and RTS2.Second, they have a fairly simple acceptance strategy that is straightforward to implement.As pointed out by Namazi et al. (2018), the authors had empirically observed that for some problems, the LAHC algorithm, unfortunately, behaves in a similar manner to the HC algorithm and does not accept worsening candidates.Thus, we utilize the DLAS algorithm to search for designs with two blocking factors.

A modification of diversified late acceptance search
Algorithm 1 provides a general scheme of a modification of the DLAS algorithm (MDLAS), which is an abridged version of Web Appendix A. Our algorithm requires four inputs: () a given design (  ), () the number of iterations ( 0 ), () the length of a fitness array ( 1 ), and () the length of a visited candidate array ( 2 ).First, in code lines 5-8, we initialize a candidate design (  ), a design (  ) that stores the best design that has been found so far, a fitness array , a visited candidate array , and the number of occurrences ( ϕ ) of ϕ max in , where ϕ max is the maximum value of .Code lines 9-13 perform  0 iterations to find an optimal design.At each iteration, the algorithm achieves three major goals.First, it finds a potential candidate design (  ), which is the best neighbor design of a given design after eliminating those that have its design efficiencies in .A neighbor design is obtained from a given design by applying a swapping mechanism to this design.Here, we consider a swapping mechanism whereby we exchange two treatments in , while fixing  (see Web Appendix B for more details).Second, this potential candidate design is accepted if the  values of it and   are identical or the  value of   < ϕ max (see Namazi et al. 2018, Section 3.1 and Web Appendix A for more details).Then, we update the temporary best design (  ) if possible and .Third, we update  using the replacement strategy (see Namazi et al. 2018, Section 3.2 and Web Appendix A for more details).We recalculate the maximum value of the array ϕ max and update the number of occurrences of ϕ max in .As shown in Algorithm 1, the major cost here is to find the best neighbor design which requires to evaluate the -efficiencies of all neighbor designs of   , denoted as  (  ).The new update formula proposed in Section 4 can be used to speed up the calculation of the -efficiencies of  (  ).But for some design parameters, the number of all neighbor designs of   () might be sizeable.To tackle this, we introduce another parameter 0 <  3 ≤ 1 to select a random sample of all neighbor designs ( ()).This allows increasing diversification and reducing the computing time when  is massive.This strategy is often used in local search algorithms.In this case, the total number of neighbor designs used our by algorithm is equal to ⌊ 3 ⌋, where ⌊⌋ is the largest integer less than or equal to .

A NEW UPDATE FORMULA
Let   and   be the treatment information matrices before and after swapping two treatments.Let Ω  and Ω  be the Moore-Penrose inverses of these matrices.Let  =   −   .This implies   =   − .Because  is a symmetric matrix, it can be factorized as Λ ⊤ , where  is the square  ×  matrix whose columns are the eigenvectors of , and Λ is the diagonal matrix whose diagonal elements are the corresponding eigenvalues.According to Henderson and Searle (1981), provided the column space of  = Λ ⊤ is a subset of that of   , or equivalently,   Ω  Λ ⊤ = Λ ⊤ , Ω  can be expressed as: where  + is the Moore-Penrose inverse of , Ψ = Λ −1 −  ⊤ Ω  , and Λ −1 is the inverse of Λ.This is a generalized version of the well-known Sherman-Morrison formula in linear algebra for the Moore-Penrose inverse of the sum of two matrices.
Since  is a binary matrix and every row has exactly one non-zero element (i.e., a treatment), by permuting rows and columns of , the general form of  can be rearranged as: where  is a  ×  permutation matrix,  1 and  2 are scalars,  is a vector of size  − 2, and  is a ( − 2) × ( − 2) zero matrix (John, 2001).Without loss of generality, it is assumed here that the first two rows in  correspond to treatments 1 and 2, and are swapped.Our challenge is to obtain the eigenvalues and the eigenvectors of the matrix   .Since the rank of this matrix is less than or equal to 2, there are at most two non-zero eigenvalues  1 and  2 , with eigenvectors  1 and  2 (John, 2001): , and If one of the two eigenvalues is equal to zero, the corresponding eigenvector is eliminated.Note that because the matrix   is a permutated version of the matrix  and is symmetric, the eigenvalues of these matrices are identical, and the eigenvectors of the matrix  are obtained by permuting the elements of the eigenvectors of the matrix   .Therefore, in order to find the eigenvalues and eigenvectors of , it is sufficient to find the eigenvalues and eigenvectors of   .
For our design problems, since  1 =  2 or  2 =  2 (see Web Appendix C for more details), we propose an alternative solution for the two eigenvectors of   : Proposition 2. The two eigenvectors of   are given by These two eigenvectors need to be normalized as , where The proof of Proposition 2 is obtained simply by showing that     =     ( = 1, 2).This proposition also holds for the case  1 ≠  2 and/or  2 ≠  2 .Note that, given an eigenvalue of a matrix, there exist many different eigenvectors before normalization (Bronson, 1991, Chap. 6).However, these eigenvectors are equivalent after normalization (possibly with different sign (±)).So, if  1 ≠  2 (or  2 ≠  2 ), the eigenvectors identified in Equations ( 11) and ( 12) are indeed equivalent after normalization.Our alternative solution for the two eigenvectors is considerably simpler than the solutions found by John (2001), requiring less calculation (i.e., there is no need to calculate  ⊤ ) and allowing the calculation of the two eigenvectors when  1 =  2 or  2 =  2 .Furthermore, the lengths of our proposed eigenvectors are easy to compute (Equation ( 13)).For these reasons, it is easy to evaluate the -efficiencies of many designs in parallel using our update formulae which are a combination of Equation ( 8), Equation (10), and Proposition 2. Using the new update formulae, the calculation of -efficiencies of many designs can be implemented in Julia using matrix computations, which leverage the advantage of Basic Linear Algebra Subprograms (BLAS) that uses both vectorization and parallelization (BLAS, 2020) (see Algorithm 2 in Web Appendix B).As illustrated in the next section, this is extremely useful when we have a large number of neighbor designs.

COMPUTATIONAL EXPERIMENTS
We first evaluate the performance of the proposed update formulae using the six examples shown in Table 1.Examples 2-6 were used in Vo-Thanh and Piepho (2020), who included a third blocking factor, in addition to row and column blocking factors, to guarantee an even distribution of cells with checks (i.e., check plots).Here, we only consider row and column blocking factors.Second, we apply our search framework to search for designs for these examples.Then, we compare our designs to those obtained by CycDe-sigN (VSNi, 2020), DiGGer (Coombes, 2020), and OPTEX (Pereira & Tobias, 2015;Piepho, 2015).Our algorithm was implemented in Julia version 1.5 under Windows 10, and our computational experiments used an Intel (R) Core (TM) i7-6600U CPU @ 2.6 GHz, and 16 GB of memory (see JCode.zip in Web Material).

Performance of the new update formulae
Given a random design () for each example in Table 1, we compare the average computing times for -efficiencies of designs in  () under four different approaches.The first one is a standard approach (  1 ) in which we compute  for each design in  () sequentially without update formulae.In this approach, we only need to recompute  in Equation ( 4) when evaluating  for each design in  () since  1 and  2 do not change.The second approach (  2 ) uses John (2001), while the third approach (  3 ) uses our proposed solution.The fourth approach (  4 ) uses our proposed solution and Algorithm 2 in Web Appendix B, which computes all neighbor designs in parallel.
The average computing times required to evaluate the -efficiencies of all designs in  () under the four approaches are shown in Table 2.As expected, the   1 approach had the longest computing time, ranging from 0.090 to 13792.997s, whereas the   2 approach took between 0.016 to 120.887 s, the   3 approach took 0.008 to 112.955 s, and the   4 approach took < 0.000 to 0.773 s.The   3 approach is slightly faster than the   2 approach.This is because the   3 approach does not require  ⊤  to be computed.The   4 approach is the fastest.Remarkably, for Example 6, this approach required 0.773 s to compute all neighbor designs, while the   3 approach required 112.955 s.This shows that there is a substantial benefit when implementing our proposed solution with parallel computing using Algorithm 2 in Web Appendix B.

Performance of the search framework
To find ARCDs with the design parameters shown in Table 1, we applied Algorithm 1 with  0 = 300,  1 = 30, TA B L E 3 Overall efficiency factor () and overall weighted efficiency factor ( Ā) of designs obtained by MDLAS, designs obtained by CycDesigN in two design menus (i.e., non-resolvable row-column designs with unequal replications and -rep single location denoted by CD1 and CD2), designs obtained by rcDiGGer, and designs obtained by OPTEX with design parameters in Table 1  Note: "-" stands for "not available".Parameters used in Algorithm 1 :  0 = 300;  1 = 30;  2 = 60;  3 = 1.MLDAS uses  as a main search criterion.
CycDesigN uses Ā as a main search criterion.rcDiGGer uses V as a main search criterion.Designs obtained by minimizing V and by maximizing  are equivalent.
OPTEX uses -optimality as a main search criterion.
A subscript asterisk (*) marks the best -efficiency of an optimal design found for each example.
and  2 = 60.For the six examples, we use  3 = 1 to test the limit of our algorithm.But in practice, we could use a smaller value of  3 to get a good solution within reasonable computing time, as will be shown later in the paper.These values ( 0 = 300,  1 = 30, and  2 = 60) were found after we performed a one-factor-at-a-time method with Algorithm 1. Table 3 shows the values of , Ā, and the average computing times required to obtain the ARCDs.
In Table 3, we also report the  and Ā values of designs obtained from CycDesigN, DiGGer, and OPTEX.For CycDesigN, we found augmented row-column designs under two design menus.Designs obtained in CycDe-sigN under the non-resolvable row-column design menu with unequal replications and under the -rep single location menu are denoted as CD1 and CD2, respectively.For Example 1, in both menus, CycDesigN obtained a design with high overall weighted efficiency factor ( Ā) with one or two tries.For Examples 2-6, we obtained designs with high overall weighted efficiency factor ( Ā) after ten tries.Often, CycDesigN converges to a good solution very fast, ranging from a few seconds to 200 s.After that, the package cannot improve the good design further.For DiGGer, we obtained designs from the command rcDiGGer in which we set high variance ratios for row and column effects, and optimized  due to the fact that DiGGer employs linear mixed models (Searle et al., 2009).Assigning high variance ratios for row and column effects should make the search work in a similar way to fixed row and column effects.rcDiGGer converges to a good solution in times ranging from a few seconds to 1278 s.After that, the package cannot improve the good design further.Detailed settings used to obtain these designs can be found in the Web Material of this paper (see JCode.zip).Because designs obtained by minimizing  and maximizing  are equivalent, we report the values of  instead of  in Table 3.For designs obtained from OPTEX, the average computing times required for the six examples range from a few seconds to 3,780 s with 100 iterations.The SAS script to obtain these designs can be found in the Web Material of this paper (see JCode.zip).Note that designs found by OPTEX are optimized according to the -optimality criterion.It is also possible to use the -optimality criterion in OPTEX to find ARCDs.However, the -optimality criterion defined in OPTEX is based on the inverse of a treatment information matrix that has a full rank, which is different from our -efficiency definition based on the Moore-Penrose inverse of the treatment information matrix that is not of full rank (see Equation ( 3)).Furthermore, the -optimality criterion in OPTEX is not invariant to nonsingular recoding of the design matrix as opposed to the -optimality criterion.For these reasons, we used the -optimality criterion in OPTEX to find ARCDs.These designs perform well in terms of both  and Ā (see in Table 3).
As shown in Table 3, given a random starting ARCD, the average computing times required to obtain a good ARCD with MDLAS range from a few seconds to 351 s.For Example 1, the designs obtained from the three approaches (i.e., CycDesigN, DiGGer, and OPTEX) and our design are equally good in terms of  and Ā.Compared to designs obtained from CycDesigN (under the two menus) for Examples 2-6, our designs are slightly better in terms of both  and Ā.For Example 6, we were able to find an augmented row-column design, while CycDesigN's non-resolvable row-column design menu with unequal replications did not find a design.For Examples 1 and 2, our designs and those obtained from DiGGer are equally good in terms of  and Ā.For Examples 3, 4, and 6, our designs are slightly better than those obtained from DiGGer.For Example 5, the design obtained from DiGGer is slightly better than ours.For Examples 1 and 2, the OPTEX designs are equally good as ours, while for Examples 3-5, our designs are slightly better than those obtained by OPTEX.As a side note, it might be possible to obtain better designs with CycDesigN, DiGGer, and OPTEX when using different starting designs.The same goes for our Algorithm 1 where we use two different random seeds when generating ARCDs.However, our limited results show that the performance of our approach is competitive compared with the existing packages in terms of both  and Ā.Furthermore, our algorithm was developed in an open-source software as opposed to the existing packages.Thus, it is accessible for those who want to use or extend our package.Finally, our algorithm has an excellent computing time and can handle various parameters for designs with two blocking factors and a complex blocking structure that cannot be handled with the existing packages as illustrated in the next section.
As mentioned before, it is possible to use a small value of  3 to get a good design within a reasonable computing time.To illustrate this point, consider Example 6, where we used  3 = 0.5.The optimal design for this example was found after 190.882 s, which is approximately 54.45% of the average computing time required with  3 = 1 (350.58s).The values of  and Ā of this design are 0.36994 and 0.41150.Compared to the optimal design with  3 = 1 ( = 0.37064 and Ā = 0.41203), we lost 0.19% for  and 0.13% for Ā.These losses are minor.However, we saved 45.55% of the computing time in getting a good design for this example.

APPLICATIONS
Here, we show four applications of our work.First, we provide an example where our proposed eigenvectors work, while the formulae in John ( 2001) cannot be used.This first application can be found in Web Appendix C. Next, we use our search framework to search for optimal designs for Examples I, II, and III.

An optimal augmented row-column design for Example 𝐈
The optimal design for Example I found by our search framework is shown in Figure 2. The average computing time required to obtain a good design this example is less than one second.The choice of the number of entries and the six replicates for the three checks enable the nice symmetry of this design.Every row (column) has three checks, denoted by 1, 2, and 3.This implies each row (column) is a complete block for the checks.Therefore, row (column) effects can be well captured via checks.

An optimal two-phase design for Example 𝐈𝐈
The average computing time required to obtain a good design for Example II is less than one second.The optimal design is shown in Figure 3.It is universally optimal under model (1): see Morgan and Uddin (1993, Theorem 1).As explained in Section 2.3, rows and columns in the design in Figure 3 correspond to field blocks and days, respectively.A cell formed by the intersection of a row and a column in this design has either zero experimental units (i.e., the cell is empty) or two experimental units (i.e., the cell is filled).Note that when removing empty cells in Figure 3, treatments are arranged in five superblocks of size 2 × 3. A row in each superblock corresponds to a field block, while a column in each superblock corresponds to a day.Both rows in a given superblock have the same set of treatments, which is one of the key features of the design.This pattern reflects the constraint that the numbers of rows and columns in every superblock divide the number of replicates, which is 6 in this example (Morgan & Uddin, 1993, Theorem 1).Note that it may be possible that for certain values of the design parameters, this constraint is violated.However, our search algorithm is still able to search for optimal two-phase designs for those settings.For example, consider the case where  = 50,  = 2,  1 = 20,  1 = 5,  2 = 25, and  2 = 4.Under these settings, the number of rows in a superblock is equal to  1 = 4, which does not divide  = 2. Using our search framework, we obtained an optimal design for this case in less than 3 s (see JCode.zip).
F I G U R E 3 An optimal two-phase design for Example II.Rows are blocks in Phase 1, while columns are blocks in Phase 2. When removing empty cells in this figure, treatments are arranged in  = 5 superblocks.Each superblock has  1 = 2 rows and  2 = 3 columns.Each intersection of a row and a column is a superplot with  = 2 experimental units.

F I G U R E 4
An optimal incomplete row-column design for Example III.An empty cell is an unusable plot.

An optimal incomplete row-column design for Example 𝐈𝐈𝐈
The optimal design for Example III is shown in Figure 4.The average computing time required to obtain a good design for this example is less than one second.Note that we would expect to have two replicates for every treatment in every row if there were no missing plots.As shown in Figure 4, in any row, most of the treatments have at most two replications.Furthermore, the -efficiency of this design is equal to 0.96498.It is very close to 0.98226, which is an upper bound for the average efficiency factor of a row-column design with no missing plots (John & Williams, 1995).This indicates that our search framework is efficient in finding a good design for an experiment where a few plots cannot be used.

CONCLUSION
In this paper, we first proposed a new search framework to obtain designs with two blocking factors (e.g., augmented row-column, two-phase, and incomplete row-column designs).Second, we proposed a new solution when the update equations of John ( 2001) cannot be used.Third, we briefly explained how to apply our proposed solution to compute the -efficiencies of many neighbor designs in parallel within a reasonable computing time when searching for optimal designs with two blocking factors.Furthermore, we found our search framework to be competitive for cases that can be tackled with existing software, while allowing more challenging problems to be dealt with.We showed that good two-phase designs and incomplete row-column designs can found under our search framework.Throughout this paper, we treated block effects as fixed and assumed var() =  2   .In agricultural applications, it is common to consider the effects of blocks in model (1) as random for recovery of inter-block information.Furthermore, experimenters often have a prior knowledge about variance components.Therefore, it would be useful to extend our proposed solution to cover this scenario.Given variance-covariance structures for the block effects and the residual error, the matrix  is known.Thus, the matrix  still has a form such as the one shown in Equation ( 9).The major challenge then is to find the three components of  (i.e.,  1 ,  2 , and ) via the matrices  and .Second, ARCDs obtained in Section 5 might not have an even distribution for check plots.Recently, Piepho et al. (2016Piepho et al. ( , 2018) ) proposed two criteria to obtain evenness of distribution of treatment replications in row-column designs, while the work of Vo-Thanh and Piepho (2020) proposed a new design class for augmented row-column experiments that can also guarantee an even distribution for check plots.Thus, it would be interesting to use our search algorithm in combination with the two criteria proposed in Piepho et al. (2016Piepho et al. ( , 2018) ) to search for ARCDs and to extend our search algorithm to the design class proposed in Vo-Thanh and Piepho (2020).Finally, it is possible to also use our algorithm with other design problems such as spatial designs and -rep designs (Cullis et al., 2006).

A C K N O W L E D G M E N T S
We would like to acknowledge the support of the Deutsche Forschungsgemeinschaft (DFG) via grant PI 377/19-2.We are further grateful to an associate editor and two anonymous referees for their invaluable comments, and Neil Coombes for his help with the DiGGer package.
Open access funding enabled and organized by Projekt DEAL.

D ATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings in this paper are available in the Supporting Information of this paper.
Nha Vo-Thanh https://orcid.org/0000-0002-9901-0545R E F E R E N C E S Williams, E.R. & Piepho, H.-P. (2015) Optimality and contrasts in block designs with unequal treatment replication.Australian & New Zealand Journal of Statistics, 57, 203-209.S U P P O R T I N G I N F O R M AT I O N Web Appendices A, B, and C referenced in Sections 3, 4, 6, as well as software codes used in this paper are available with this paper at the Biometrics website on Wiley Online Library.How to cite this article: Vo-Thanh, N. & Piepho, H.-P. (2023) Generating designs for comparative experiments with two blocking factors.Biometrics, 79, 3574-3585.https://doi.org/10.1111/biom.13913 along with the average computing time in seconds required by MDLAS for one random starting design.