Design of Order-of-Addition Experiments

Design of Order-of-Addition Experiments Jiayu Peng * 1 , Rahul Mukerjee †2 , and Dennis K.J. Lin ‡1 1 Department 2 Indian of Statistics, Penn State University Institute of Management Calcutta Abstract In an order-of-addition experiment, each treatment is a permutation of m components. It is often unaffordable to test all the m! treatments, and the design problem arises. We consider the model in which the response of a treatment depends on the pairwise orders of the components. The optimal design theory under this model is established, and the optimal values of the D-,A-, E-, and M.S.-criteria are derived. We identify a special constraint on the correlation structure of such designs. The closed-form construction of a class of optimal designs is obtained, with examples for illustration. Key Words: Design equivalence; Optimal design; Pairwise order. 1 Introduction In many chemical experiments, a number of reactants are added into the system sequentially rather than simultaneously. The formation, e.g., the amount/size/purity of the reaction product often depends on the order-of-addition (OofA) of different reactants. Fuleki and Francis (1968) is one of the early work that addresses the OofA issue: “The order of addition of the lead acetate (before or after the pH adjustment) had a definite influence on the reaction. Higher recoveries were obtained by adjusting the pH after lead acetate addition.” During the past decades, the OofA effect is frequently mentioned in the area of bio-chemistry (Shinohara and Ogawa 1998), food science (Jourdain et al. 2009), nutritional science (Karim et al. 2000), and pharmaceutical science (Rajaonarivony et al. 1993), etc. The general consideration of OofA is about the addition of any m components, not necessarily chemical reagents. Ultimately, an OofA experiment is attempt to find the optimal addition order. As an example in genetics, when people construct the phylogenetic trees (a.k.a. evolution trees), the property of the obtained tree depends on the order-of-addition of taxa (Olsen et al. 1994). Since an exhaustive search of all sequences is not affordable, often * [email protected] † ‡ [email protected] [email protected] 1 a number of randomly selected taxa orders are tested (Stewart et al. 2001). However, randomly selected orders may not be informative; rather, it is desired to carefully design a subset of orders to be tested. Van Nostrand (1995) appears to be the first statistical reference on the OofA design. He suggested to consider the design whose goal is to detect the pairwise-order (PWO) effects. In fact, the pairwise-order effects are of great interest in many experiments: practitioners often expect conclusions like “Adding A before or after B has a significant influence on the response” (see the aforementioned example in Fuleki and Francis 1968). Based on Van Nostrand’s PWO model, Voelkel (2017) proposed a number of design criteria. He found some efficient PWO designs of which the correlation structures (moment matrices) are the same as the full PWO designs. These designs are intuitively optimal, e.g., under the D-optimality criterion, but there is no theoretical support. In this work, we show that the moment matrix of the full PWO design is D-, A-, E-, M.S.-, and G-optimal. Moreover, an arbitrary PWO design is D-, A-, or M.S.-optimal if and only if it has the same moment matrix as the full PWO design. We prove the optimality results via two different approaches, both of theoretical interest. The first one, a non-computational proof, is based on the approximate theory (Kiefer 1974), making use of the concavity and signed permutation invariance of the optimality criteria. The second one, a computational approach, is via identifying a constraint on the correlations of all PWO designs, then identifying a constraint on the eigenstructure of their moment matrices. The remainder of this paper is as follows. Section 2 formulates the model and defines the PWO design. Section 3 proves the optimality results on PWO designs, based on the approximate theory. In Section 4, we evaluate the moment matrix of the full PWO design, and derive the exact values of the benchmarks for assessing any other PWO design. Section 5 presents a special constraint on the correlation structure of any PWO design. This constraint, derived from the transitive property, intuitively implies how orthogonal a PWO design could be. Based on this constraint, we briefly introduce our second approach to prove the optimality results stated in Section 3. Section 6 is devoted to a closed-form construction for a class of optimal fractional PWO designs. The result confirms that for any number of components (m), there exist fractional PWO designs (with number of runs as a fraction of m!) that attain the optimal moment matrix. Section 7 discusses the potential applications in constructing efficient PWO designs with smaller number of runs. For simplicity of the presentation, all proofs are deferred to the Appendix. 2 Model Formulation Suppose there are m (≥ 3) components 1, 2, . . . , m, which can be ordered in m! ways. Any such ordering (say a = a1 a2 ⋯am ), which is a permutation of 1, . . . , m, is a treatment. Denote the set of all the m! treatments by A. For a ∈ A, write τ (a) for the treatment effect of a, i.e., τ (a) is the expectation of any observation arising from treatment a. As usual, it is assumed that the observations have equal variance and are uncorrelated. Under the first-order model, for every 2 a ∈ A, we have τ (a) = β0 + ∑ zjk (a)βjk , (1) jk∈S where S is set of all pairs jk for 1 ≤ j < k ≤ m, βjk ’s and β0 are unknown parameters, and for each jk ∈ S, ⎧ ⎪ if j precedes k in a, ⎪1 zjk (a) = ⎨ (2) ⎪ −1 if k precedes j in a. ⎪ ⎩ As an example, when m = 4, S = {12, 13, 14, 23, 24, 34}, and for the treatment a = 2143, z12 (a) = −1, z13 (a) = +1, . . . , and z34 (a) = −1. Let β̃ = (β12 , β13 , . . . , β(m−1)m )T , where T denotes the T transpose. Then β = (β0 , β̃ T ) represents the parametric vector of interest. Similarly, for any T T a ∈ A, let z(a) = (z12 (a), z13 (a), . . . , z(m−1)m (a)) and x(a) = (1, z(a)T ) . Then (1) can be expressed as τ (a) = β0 + z(a)T β̃ = x(a)T β (3) Write q = (m2 ) and p = q + 1. Then β̃ and z(a) are q × 1 while β and x(a) are p × 1. The design of OofA experiments is to choose a collection D = {a1 , a2 , . . . , aN }, where ai ∈ A for each i. We represent each D by a pairwise-order (PWO) design, defined as Z = [z(a1 ), z(a2 ), T vf b . . . , z(aN )] . Then Z is the design matrix in model (1). Let Df be the full design which replicates each treatment in A once, and let Zf be the full PWO design, i.e., the PWO design of Df . For example, below are the Df and Zf for m = 3. Each row in Table 1 (a) represents a treatment in Df . Table 1: The full OofA design and full PWO design with 3 components (b) Full PWO design Zf (a) Full OofA design Df 1 1 2 2 3 3 2 3 1 3 1 2 z12 z13 z23 + + + + + − − + + − − + + − − − − − 3 2 3 1 2 1 Remark 1. Due to the transitive property of order, the region of PWO design points does not include all the level combinations. For example, when m = 3, the points (+, −, +) and (−, +, −) are not valid for a PWO design. 3 Theory of Optimal PWO Design Given any m ≥ 3, consider the full OofA design Df , which replicates each treatment once. Let Zf be the full PWO design. For inference on β, we first explore the optimality of Zf under commonly 3 used criteria among all designs with the same number, m!, of runs. Although Zf is impractical for larger m, any optimality result will be useful as a benchmark for assessing smaller designs. Such optimality is anticipated because, with run size m!, a design that replicates some treatments more than once while omitting some others altogether is intuitively unappealing. However, unlike in other similar situations like traditional full factorials, the information matrix of Zf , obtained later in Section 4, is complicated. For example, it does not possess complete symmetry (compared with Kiefer 1975) in the sense of having all off-diagonal elements equal. As such, a direct combinatorial proof of the optimality of Zf is quite challenging. An alternative approach, based on the approximate theory, yields a subtle non-computational proof that does not require explicit evaluation of the information matrix of Zf . 3.1 Approximate Theory Background and Optimality Criteria Consider an N −run OofA design D, where any treatment a ∈ A is replicated r(a) (≥ 0) times, with the integers r(a) summing to N . By (3), the PWO design of D (i.e., Z) has information matrix M(w) = ∑ w(a)x(a)x(a)T , (4) a∈A where w(a) = r(a)/N . In approximate theory, the requirement on the design weights w(a) to be integral multiples of 1/N is relaxed, i.e., each w(a) (a ∈ A) is allowed to be any non-negative quantity, subject to ∑a∈A w(a) = 1. Then w = {w(a) ∶ a ∈ A} is called a design measure, having moment matrix M(w) as shown in (4). In particular, the full design Zf corresponds to the uniform design measure w0 over A which has moment matrix Mf = M(w0 ) = 1 ∑ x(a)x(a)T . m! a∈A (5) Let M denote the class of p × p non-negative definite matrices. Recall that a signed permutation matrix is a square matrix having exactly one non-zero entry in each row and column, where any non-zero entry is either 1 or −1. We consider optimality criteria φ(⋅) that are (i) concave over M and (ii) signed permutation invariant, i.e., φ(RT MR) = φ(M) for every signed permutation matrix R of order p and every M ∈ M . Given any such criterion φ(⋅), a design measure w is called φ-optimal if it maximizes φ(M(w)) among all design measures. The commonly used D, A-, and E-criteria correspond to φ(M) = log det(M), − tr(M−1 ), and λmin (M) respectively, where log det(M) and − tr (M−1 ) are interpreted as −∞ for singular M, and λmin stands for the smallest eigenvalue. Because all design measures here have the same tr (M(w)) (= p) by (4), our framework also covers the M.S.-optimality criterion of Eccleston and Hedayat (1974), with φ(M) = − tr(M2 ). This is equivalent to the E(s2 )-criterion for two-level supersaturated designs (Booth and Cox 1962). All these criteria are concave and signed permutation invariant. 3.2 Optimality of the Uniform Design Measure We now establish the optimality of the uniform design measure w0 , or equivalently, that of the full PWO design Zf . Our main idea is as follows. Although the space of PWO design points seems 4 irregular (see Remark 1), it is symmetric in the following sense: there exists a group action on the design space, which is essentially induced by the symmetric group formed by all the treatments in A. Technically, the action of each element can be represented as a signed permutation matrix (see Lemma 1 in the appendix). Then making use of the signed permutation invariance as well as the concavity of optimality criteria, we derive the following result. Theorem 1. The uniform design measure w0 is φ-optimal for every optimality criterion φ(⋅) which is concave and signed permutation invariant. Theorem 1 has several important implications: (i) It shows, in particular, the D-, A-, E-, and M.S.-optimality of the uniform design measure w0 . In view of the equivalence theorem, the D-optimality of w0 also implies its G-optimality, i.e., max x(a)T M−1 (w)x(a) ≥ p = max x(a)T M−1 f x(a) a∈A (6) a∈A for every design measure w; see, e.g., Silvey (2013), Chapter 3 for more details. By (6), w0 minimizes the maximum variance of the estimated responses at a ∈ A. (ii) The values of these criteria for the full design Zf (as detailed in the next section) provide benchmarks for evaluating any (not full) PWO design. (iii) While w0 may not be the unique design measure enjoying the above optimality properties, any other design measure which is D-, A-, or M.S.-optimal must have the same moment matrix Mf as w0 . This is because the φ(⋅) functions, corresponding to these three criteria, are strictly concave on the manifold M . (iv) A design with smaller number of runs is D-, A-, E-, M.S.-, and G-optimal if its moment matrix is the same as Mf . Designs of this kind are of special interest to us and will be discussed in Section 6. 4 Properties of the Full PWO Design Even though the optimality proof indicated in Section 3 does not require explicit knowledge of Mf , we will need to find Mf and its eigenvalues for comparative purposes, i.e., to assess the efficiency of a given design measure or given exact design, under various optimality criteria. The form of Mf can be presented as below. ̃f is ̃f ), where M Proposition 1. The full PWO design Zf has the moment matrix Mf = diag(1, M a q × q matrix with rows and columns indexed by the elements of S, such that for any j1 k1 , j2 k2 ∈ S, ̃f is the (j1 k1 , j2 k2 )-th element of M ⎧ 1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪1/3 ̃f (j1 k1 , j2 k2 ) = ⎪ M ⎨ ⎪ −1/3 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩0 if j1 = j2 and k1 = k2 if j1 = j2 , k1 ≠ k2 or j1 ≠ j2 , k1 = k2 if j1 = k2 or k1 = j2 otherwise. 5 (7) For instance, if m = 4, then S = {12, 13, 14, 23, 24, 34}, and ⎡ 1 1/3 1/3 −1/3 ⎢ ⎢ 1/3 1 1/3 1/3 ⎢ ⎢ ⎢ 1/3 1/3 1 0 ̃f = ⎢ M ⎢ −1/3 1/3 0 1 ⎢ ⎢ ⎢ −1/3 0 1/3 1/3 ⎢ ⎢ 0 −1/3 1/3 −1/3 ⎣ −1/3 0 1/3 1/3 1 1/3 0 −1/3 1/3 −1/3 1/3 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥, ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ̃f (12, 13) = 1/3, M ̃f (12, 14) = 1/3, M ̃f (12, 23) = −1/3, and so on. If ̃f (12, 12) = 1, M where M we interpret any zjk (a) as the level of a factor Zjk in the treatment a, then Proposition 1 implies that: (i) A full PWO design is always balanced, that is, each factor is run the same number of times at the + and − levels. (ii) In the full design, the correlation is 1/3 between factors Zjk and Zjl for any j and k ≠ l, both of which indicating whether j precedes the other component - this is called a synergistic pair of PWO factors in Voelkel (2017). The correlation is −1/3 between factors Zkj and Zjl for any k < j < l - this is called an antagonistic pair of PWO factors. If two PWO factors do not involve any common component, they are uncorrelated in the full design. Based on Proposition 1, we can derive the eigenvalues of Mf and hence the explicit values of the various design criteria for the full design. ), Theorem 2. Mf has eigenvalues 1, (m + 1)/3 and 1/3, with multiplicities 1, m − 1 and (m−1 2 respectively. Then for the full design: (m + 1)m−1 p [det(Mf )] = [ ] , 3q 3m(m − 1)2 ) , tr (M−1 = 1 + f 2(m + 1) 1 λmin (Mf ) = , and 3 (m − 1)m(2m + 5) . tr (M2f ) = 1 + 18 1 1/p ̃f . Consider the q−dimensional space indexed by Next, we describe the eigenvectors of M the elements in S. Let ejk be the vector of which the jk-element is 1 and other elements are 0. For example, if m = 4 then S = {12, 13, 14, 23, 24, 34}, and e12 = (1, 0, 0, 0, 0, 0)T , . . . , e34 = (0, 0, 0, 0, 0, 1)T . For any 1 ≤ l ≤ m and any jk ∈ S, let δl,jk ⎧ 1 ⎪ ⎪ ⎪ ⎪ = ⎨−1 ⎪ ⎪ ⎪ ⎪ ⎩0 if l = j if l = k otherwise, 6 (8) and for any 1 ≤ l ≤ m, define ul = ∑ δl,jk ejk . (9) vjkl = ejk − ejl + ekl , (10) jk∈S Further define T = {jkl ∶ 1 ≤ j < k < l ≤ m}, and let for any jkl ∈ T . For example, when m = 4, u1 = e12 + e13 + e14 , u2 = −e12 + e23 + e24 , u3 = −e13 − e23 + e34 , and u4 = −e14 − e24 − e34 . On the other hand, v123 = e12 − e13 + e23 , v124 = e12 − e14 + e24 , v134 = e13 − e14 + e34 , and v234 = e23 − e24 + e34 . We have: Proposition 2. ̃ f = m + 1 PU + 1 PV , M 3 3 where U = span{ul ∶ 1 ≤ l ≤ m}, V = span{vjkl ∶ jkl ∈ T }, and for any linear space W , PW indicates the projection matrix onto W . 5 The Constraint on the Correlation Structure of PWO Designs Although not required for proving the optimality result (Section 3), some knowledge about the explicit form of M for any (not full) design helps us in intuitively understanding how orthogonal a PWO design could be. Consider any N -run, m-component OofA design D and its moment matrix ̃ = 1 ∑a∈D z(a)z(a)T , i.e., M ̃ is obtained by removing the first row and column of M. M. Let M N ̃ is indexed by the elements in S, and M(j ̃ 1 k1 , j2 k2 ) indicates the correlation between the Then M ̃ has the following property. PWO factors Zj1 k1 and Zj2 k2 . M ̃ For any 1 ≤ j < k < l ≤ m, Proposition 3. Consider any PWO design Z and its corresponding M. we have: ̃ ̃ ̃ M(jk, jl) − M(jk, kl) + M(jl, kl) = 1. (11) Note that Proposition 3 is derived from the transitive property of order, which says: For any three different components j1 , j2 , and j3 , if j1 precedes j2 and j2 precedes j3 , then j1 must precede ̃ ̃ ̃ j3 . From Proposition 3 we get M(jk, jl)2 + M(jk, kl)2 + M(jl, kl)2 ≥ 1/3. The equality holds if ̃ ̃ ̃ ̃=M ̃f . Thus M ̃f and only if M(jk, jl) = −M(jk, kl) = M(jl, kl) = 1/3, which holds when M ̃ is intuitively the most efficient among all M’s. Also, Proposition 3 implies that the full design ̃ see, e.g., minimizes the rmax , i.e., the maximum absolute value of the off-diagonal elements of M; Liu and Dean (2004). An interesting result following Proposition 3 is stated below. ̃f with eigenvalue 1/3 (see Section 4). Corollary 1. Let vjkl ’s (jkl ∈ T ) be the eigenvectors of M ̃ Then for any jkl ∈ T and the M corresponding to any PWO design, it holds that T ̃ vjkl Mvjkl = 1. 7 (12) ̃f vjkl = (1/3)vjkl and v T vjkl = 3, equation (12) apparently holds for M ̃f . In Remark 2. Since M jkl ̃f is replaced by an arbitrary M, ̃ vjkl is (generally) no longer its eigenvector, but other words, as M ̃ stays the same. This is an essential constraint on the spectrum the norm of vjkl , with respect to M, ̃ of M. With this result, we are able to prove Theorem 1 via a computational approach (see the Appendix). 6 Closed-form Construction of a Class of Optimal Fractional PWO Designs From Section 3, it is evident that a PWO design is optimal under each of the D-, A-, E-, G-, and M.S.-criteria if it has the same moment matrix as the full design. Does such optimal design exist among designs with less than m! runs? Voelkel (2017) found several such designs for m = 4, 5 or 6, via computer search. In this section, we systematically construct a class of optimal, fractional PWO designs for any m ≥ 4. Example 1. When m = 4, the following half-fractional OofA design has the same moment matrix as the full design and hence is optimal under all the criteria mentioned above. Table 2: Optimal half-fractional OofA design with 4 components 1 2 4 3 1 3 4 2 1 4 3 2 2 1 3 4 3 1 2 4 4 1 2 3 3 4 1 2 2 4 1 3 2 3 1 4 4 3 2 1 4 2 3 1 3 2 4 1 The above design can be obtained from the Design 1 in Table 2 of Voelkel (2017), via relabeling the components and row permutation. 1 To clearly show the structure of this design, define B1 = [ 2 2 4 1 4 2 3 B2 = [ ], B3 = [ ], and B3 = [ ]. For any matrix A 4 2 4 1 3 2 3 4 1 3 2 ], B1 = [ ], B2 = [ ], 4 3 3 1 1 = (aij )1≤i≤n,1≤j≤k , define ⊖A ∶= 4 3 (ai,k+1−j )1≤i≤n,1≤j≤k . The operator ⊖ is to reverse each row of a matrix, for example, ⊖B1 = [ ]. 3 4 8 In this way we can represent the design in Table 2 as: ⎡ ⎢ B1 ⎢ ⎢ ⎢⊖B1 ⎢ ⎢ ⎢ B2 D = ⎢⎢ ⎢⊖B2 ⎢ ⎢ ⎢ B3 ⎢ ⎢ ⎢⊖B3 ⎣ ⎤ B1 ⎥ ⎥ ⎥ B1 ⎥⎥ ⎥ B2 ⎥⎥ ⎥. B2 ⎥⎥ ⎥ B3 ⎥⎥ ⎥ B3 ⎥⎦ It will now be shown how the design structure, as we observe above, entails an extension allowing us to cover the case of general even m. To rigorously describe the approach, we bring in several notations. Definition 1 (Combination matrix). Given any 2 ≤ r < m, the matrix C(m,r) denotes the collection of all r−combinations of the set {1, 2, . . . , m}. Namely, C(m,r) is an (mr) × r matrix, with each row being a subset of {1, 2, . . . , m}, and all rows being distinct. The elements in each row of C(m,r) are placed in ascending order. Different rows of C(m,r) are placed in the lexicographic order. Definition 2 (Complement of C(m,r) ). The matrix C (m,r) is an (mr) × r matrix determined by the following two conditions. First, in the columnwise-combined matrix [C(m,r) C contains all elements in {1, 2, . . . , m}. Second, the numbers in each row of C ascending order. ⎡1 ⎢ ⎢1 ⎢ ⎢ ⎢1 As an example, C(4,2) = ⎢⎢ ⎢2 ⎢ ⎢2 ⎢ ⎢3 ⎣ ⎡3 2⎤⎥ ⎢ ⎢2 3⎥⎥ ⎢ ⎥ ⎢ ⎥ ⎢2 (4,2) 4⎥ ⎢ and C = ⎢1 ⎥ 3⎥ ⎢ ⎢ ⎥ ⎢1 4⎥ ⎥ ⎢ ⎢1 4⎥⎦ ⎣ (m,r) (m,r) ], each row are placed in 4⎤⎥ 4⎥⎥ ⎥ 3⎥⎥ . 4⎥⎥ ⎥ 3⎥ ⎥ 2⎥⎦ Definition 3 (Full permutation matrix). (i) For any set {a1 , a2 , . . . , ak }, its full permutation matrix is a k! × k matrix, with each row being a permutation of a1 , a2 , . . . , ak , and all rows being distinct. (m,r) denotes the full permutation matrix (ii) Given any 2 ≤ r < m and 1 ≤ i ≤ (mr), the matrix Bi (m,r) of the ith row of the combination matrix C(m,r) . The matrix Bi (m,r) matrix of the ith row of C . (4,2) For example, B1 denotes the full permutation (4,2) (4,2) 3 4 1 2 3 4 1 2 (4,2) ], and B1 ], . . . , B6 =[ =[ ], . . . , B6 ]. =[ =[ 4 3 2 1 4 3 2 1 9 For any r ≥ 2 and m = 2r, we define ⎡ (m,r) ⎢ B1 ⎢ ⎢ (m,r) ⎢⊖B1 ⎢ ⎢ (m,r) ⎢B ⎢ 2 ⎢ (m,r) D ∶= ⎢⎢⊖B2 ⎢ ⎢ ⋮ ⎢ ⎢ ⎢ (m,r) ⎢ BL ⎢ (m,r) ⎢ ⎢⊖BL ⎣ (m,r) ⎤ B1 ⎥ ⎥ (m,r) ⎥ B1 ⎥⎥ (m,r) ⎥ B2 ⎥⎥ ⎥ (m,r) B2 ⎥⎥ , ⎥ ⋮ ⎥⎥ ⎥ (m,r) ⎥ BL ⎥ ⎥ (m,r) ⎥ BL ⎥⎦ (13) where L = (mr)/ 2 and ⊖ is the operator defined in Example 1. Treating each row in D as a treatment in the OofA design, D corresponds to an optimal PWO design with m!/r! runs. From ̃ ≥ 2r + 1 such optimal designs with m = 2r, one can easily obtain an optimal design with any m ̃ components and m!/r! runs. In summary, we have the following theorem. Theorem 3. For any m ≥ 4 and any 2 ≤ r ≤ m/2, there exist an optimal PWO design with m components and m!/r! runs. 7 Discussion We have established the optimality theory of PWO design and have made a first step in the systematic construction of optimal fractional PWO designs. Theorems 1 and 2 provide a useful benchmark allowing the search for highly-efficient PWO designs with smaller number of runs. Algorithmic methods as well as possible refinements of Theorem 3 appear to be very promising in this regard. This is currently under investigation. Our work also has potential applications as follows: 1. Theorem 1 and Proposition 1 imply that a PWO design has the same moment matrix as the full design and hence is optimal if and only if it is balanced and E(s2 )-optimal. Then in view of Xu (2003), such optimality for PWO designs is equivalent to the minimum-moment-aberration criterion (up to order 2). Consequently, algorithms based on minimum-moment-aberration such as Xu (2002) and Lekivetz et al. (2015) can be potentially employed to find efficient PWO designs. 2. Using the rowwise-reversion technique (in Section 6) as well as the columnwise-reversion, we have found several small-run optimal designs, for example, a design with m = 8 that requires only 168 runs. A generalized method is under development. 3. The optimal designs obtained in Section 6 are naturally blocked; see equation (13). One can choose a subset of such designs by selecting only one run from each block. The obtained small-run designs are generally efficient, as shown by our initial results. Future work can be done to optimize designs of such type. 10 References T. Fuleki and F. J. Francis. Quantitative methods for anthocyanins. Journal of Food Science, 33 (3):266–274, 1968. A. Shinohara and T. Ogawa. Stimulation by rad52 of yeast rad51-mediated recombination. Nature, 391(6665):404–407, 1998. L. S. Jourdain, C. Schmitt, M. E. Leser, B. S. Murray, and E. Dickinson. Mixed layers of sodium caseinate + dextran sulfate: influence of order of addition to oil-water interface. Langmuir, 25 (17):10026–10037, 2009. M. Karim, K. McCormick, and C. T. Kappagoda. Effects of cocoa extracts on endotheliumdependent relaxation. The Journal of Nutrition, 130(8):2105S–2108S, 2000. M. Rajaonarivony, C. Vauthier, G. Couarraze, F. Puisieux, and P. Couvreur. Development of a new drug carrier made from alginate. Journal of Pharmaceutical Sciences, 82(9):912–917, 1993. G. J. Olsen, H. Matsuda, R. Hagstrom, and R. Overbeek. Fastdnaml: a tool for construction of phylogenetic trees of dna sequences using maximum likelihood. Computer Applications in the Biosciences: CABIOS, 10(1):41–48, 1994. C. A. Stewart, D. Hart, D. K. Berry, G. J. Olsen, E. A. Wernert, and W. Fischer. Parallel implementation and performance of fastdnaml-a program for maximum likelihood phylogenetic inference. In Supercomputing, ACM/IEEE 2001 Conference, pages 32–32. IEEE, 2001. R. C. Van Nostrand. Design of experiments where the order of addition is important. In ASA Proceedings of the Section on Physical and Engineering Sciences, pages 155–160, Alexandria, Virginia, 1995. American Statistical Association. J. G. Voelkel. The design of order-of-addition experiments. arXiv preprint arXiv:1701.02786, 2017. J. Kiefer. General equivalence theory for optimum designs (approximate theory). The Annals of Statistics, 2(5):849–879, 1974. J. Kiefer. Construction and optimality of generalized youden designs. In J. N. Srivastava, editor, A Survey of Statistical Design and Linear Models. North Holland, Amsterdam, 1975. J. A. Eccleston and A. Hedayat. On the theory of connected designs: characterization and optimality. The Annals of Statistics, 2(6):1238–1255, 1974. K. H.V. Booth and D. R. Cox. Some systematic supersaturated designs. Technometrics, 4(4): 489–495, 1962. S. Silvey. Optimal design: an introduction to the theory for parameter estimation, volume 1. Springer Science & Business Media, 2013. 11 Y. Liu and A. Dean. k−circulant supersaturated designs. Technometrics, 46(1):32–43, 2004. H. Xu. Minimum moment aberration for nonregular designs and supersaturated designs. Statistica Sinica, 13(3):691–708, 2003. H. Xu. An algorithm for constructing orthogonal and nearly-orthogonal arrays with mixed levels and small runs. Technometrics, 44(4):356–368, 2002. R. Lekivetz, R. Sitter, D. Bingham, M.S. Hamada, L.M. Moore, and J.R. Wendelberger. On algorithms for obtaining orthogonal and near-orthogonal arrays for main-effects screening. Journal of Quality Technology, 47(1):2–13, 2015. Appendix: Proofs Let π = π1 π2 ⋯πm be an order in A, i.e., π1 π2 ⋯πm is a permutation of 1, 2, . . . , m. Note that π induces a mapping on A, namely, π(a) = πa1 πa2 ⋯πam , (14) for any a = a1 a2 ⋯am ∈ A. We write π(a) as πa for simplicity. Lemma 1 below helps in the proof of Theorem 1. Lemma 1. Given any π ∈ A, there exists a signed permutation matrix R(π) of order p, such that x(πa)T = x(a)T R(π), for every a ∈ A. Proof of Lemma 1. Given any π = π1 ⋯πm , for any jk ∈ S, let π̄j π̄k equal πj πk if πj < πk , and equal πk πj if πj > πk . Then by (2) and (14), for every jk ∈ S and every a ∈ A, ⎧ ⎪ ⎪zjk (a) zπ̄j π̄k (πa) = ⎨ ⎪ ⎪ ⎩−zjk (a) if πj < πk , if πj < πk , ̃ because j precedes k in a if and only if πj precedes πk in πa. Thus, z(πa)T = z(a)T R(π), for ̃ every a ∈ A. Here R(π) is a signed permutation matrix of order q such that, for each jk ∈ S, the ̃ (jk, π̄j π̄k )-th element of R(π) is 1 or −1, according as πj < πk or πj > πk , respectively, and all ̃ other elements of R(π) are zeros. It now follows that x(πa)T = x(a)T R(π), for every a ∈ A, ̃ where R(π) = diag(1, R(π)) is a signed permutation matrix of order p. Proof of Theorem 1. Consider a design measure w = {w(a) ∶ a ∈ A}. For any π ∈ A, let πw be the design measure that assigns, for each a ∈ A, weight w(a) on treatment πa. Because φ(⋅) is concave, 1 1 (15) φ( ∑ M(πw)) ≥ ∑ φ(M(πw)). m! π∈A m! π∈A 12 Now, for every fixed a ∈ A, note that {πa ∶ π ∈ A} = A, so that by (5), ∑π∈A x(πa)x(πa)T = ∑a∈A x(a)x(a)T = m!Mf . Hence by (4), 1 1 ∑ M(πw) = ∑ ∑ w(a)x(πa)x(πa)T = ∑ w(a)Mf = Mf . m! π∈A m! π∈A a∈A a∈A (16) Also, by (4) and Lemma 1, for any π ∈ A, there exists a signed permutation matrix R(π) such that M(πw) = ∑ w(a)x(πa)x(πa)T a∈A = R(π)T { ∑ w(a)x(a)x(a)T } R(π) = R(π)T M(w)R(π), a∈A and hence φ(M(πw)) = φ(M(w)), because φ(⋅) is signed permutation invariant. Consequently, in view of (16), from (15) we get φ(Mf ) ≥ φ(M(w)), and the result follows. ̃f eαβ = [M(jk, ̃ Proof of Theorem 2 and Proposition 2. For any αβ ∈ S, M αβ)]jk∈S . Then for ̃ ̃ f ul = M ̃f ∑jk∈S δl,jk ejk = [gαβ ] any 1 ≤ l ≤ m, M αβ∈S , where gαβ = ∑jk∈S δl,jk M(jk, αβ). By ̃ ̃ αβ). To show that αβ) − ∑j<l M(jl, the definition of δl,jk in equation (8), gαβ = ∑k>l M(lk, ̃f ul = (m + 1)/3 ⋅ ul , it suffices to show that gαβ = (m + 1)/3 ⋅ δl,αβ for each αβ ∈ S. M ̃f (lk, αβ) = 0 for any k > l (by Proposition 1). On the other hand, (1) If α < β < l, then M ̃f (jl, αβ) = M ̃f (αl, αβ) + M ̃f (βl, αβ) = 1/3 + (−1/3) = 0. Therefore gαβ = 0. ∑j<l M (2) For l < α < β, similarly we can show that gαβ = 0. ̃ ̃f (lβ, αβ) = ̃f (jl, αβ) = M ̃f (αl, αβ) = 1/3, and ∑k>l M(lk, αβ) = M (3) If α < l < β, then ∑j<l M 1/3. Therefore gαβ = 0. ̃f (jl, αβ) equals 1 if j = α, and equals 1/3 otherwise. Thus (4) If α < β = l, then M ̃f (jl, αβ) = 1 + (l − 2)/3 = (l + 1)/3. On the other hand, M(lk, ̃ αβ) = −1/3 for any k > l, ∑j<l M ̃ so ∑k>l M(lk, αβ) = −(m − l)/3. Hence gαβ = −(m − l)/3 − (l + 1)/3 = −(m + 1)/3, which equals (m + 1)/3 ⋅ δl,αβ when α < β = l. (5) For l = α < β, similarly we can show that gαβ = (m + 1)/3 ⋅ δl,αβ . ̃f with eigenvalue (m + 1)/3. In summary, every ul is a eigenvector of M ̃f vjkl = M ̃f ejk − Next, consider vjkl for any jkl ∈ T = {(j, k, l) ∶ 1 ≤ j < k < l ≤ m}. M ̃f ejl + M ̃f ekl = [hαβ ] ̃ ̃ ̃ M αβ∈S , where hαβ = M(αβ, jk) − M(αβ, jl) + M(αβ, kl). ̃f (αβ, jk) = M ̃f (αβ, jl) = M ̃f (αβ, kl) = 0. Hence hαβ = 0. (1) If {α, β}∩{j, k, l} = ∅, then M ̃ ̃ ̃f (αβ, kl). Then (2) If α = j and β ∉ {k, l}, then Mf (αβ, jk) = Mf (αβ, jl) = 1/3 and M hαβ = 0. ̃f (αβ, jk) = 1, M ̃f (αβ, jl) = 1/3, and M ̃f (αβ, kl) = −1/3. (3) If α = j and β = k, then M Then hαβ = 1/3. ̃f (αβ, jk) = 1/3, M ̃f (αβ, jl) = 1, and M ̃f (αβ, kl) = 1/3. Then (4) If α = j and β = l, M hi = −1/3. ̃f (αβ, jk) = −1/3, M ̃f (αβ, jl) = 0, and M ̃f (αβ, kl) = 1/3. Then (5) If α = k and β ≠ l, M hαβ = 0. 13 ̃f (αβ, jk) = −1/3, M ̃f (αβ, jl) = 1/3, and M ̃f (αβ, kl) = 1. Then (6) If α = k and βi = l, M 1 hαβ = 3 . ̃f (αβ, jk) = 0 and M ̃f (αβ, jl) = M ̃f (αβ, kl) = −1/3. Then hαβ = 0. (7) If α = l, M ̃f (αβ, jk) = M ̃f (αβ, jl) = −1/3 and M ̃f (αβ, kl) = 0. Then hαβ == 0. (8) If β = j, M ̃f (αβ, jk) = 1/3, M ̃f (αβ, jl) = 0, and M ̃f (αβ, kl) = −1/3. Then (9) If β = k and α ≠ j, M hαβ = 0. ̃f (αβ, jk) = −1/3 and M ̃f (αβ, jl) = M ̃f (αβ, kl) = 1/3. Then (10) If β = l and α ∉ {j, k}, M hαβ = 0. In summary, hαβ equals 1/3 if αβ = jk or kl, −1/3 if αβ = jl, and 0 otherwise. Hence ̃f vjkl = (1/3)vjkl , for any jkl ∈ T . M Finally, consider the dimensions of the two eigenspaces, U = span{ul }1≤l≤m and V = span{vjkl }jkl∈T . m We first show that dim(U ) = m − 1. Note that ∑m l=1 ul = ∑jk∈S ∑l=1 δl,jk ejk . For any jk ∈ S, m m ∑l=1 δl,jk = δj,jk + δk,jk = 0, and therefore ∑l=1 ul = 0, which implies that dim(U ) < m. On the other hand, we claim that u1 , u2 , . . . , um−1 form a space of rank m − 1. If not, there exist θ1 , θ2 , . . . , θm−1 ∈ R which are not all zero, such that θ1 u1 +θ2 u2 +⋅ ⋅ ⋅+θm−1 um−1 = 0. Now consider vector eαm for any 1 ≤ α ≤ m − 1. Note that ejk ’s are orthogonal to each other, it is easy to see that eTαm uα = 1 and eTα<m ul = 0 for any l ≠ α, 1 ≤ l ≤ m−1. Thus 0 = eTαm (θ1 u1 + ⋅ ⋅ ⋅ + θm−1 um−1 ) = θα , for any 1 ≤ α ≤ m − 1, contradiction! Thus u1 , u2 , . . . , um−1 are linearly independent, and then dim(U ) = m − 1. Next consider dim(V ). For each 3 ≤ t ≤ m, define Wt = span {vjkl ∶ l = t, jkl ∈ T }. Namely, W3 = span{e12 − e13 + e23 }, W4 = span {e12 − e14 + e24 , e13 − e14 + e34 , e23 − e24 + e34 }, so on so ̃t ∶= Wt − ⋃t−1 ̃ ̃ ̃ ̃ forth. Denote W s=3 Ws (and W3 = W3 ), so that V is a disjoint union of W3 , W4 , . . . , Wm . ̃3 ) = 1. We now show that dim(W ̃t ) = t − 2 for any t ≥ 4. Obviously, dim(W (1) Consider the following elements in Wt : wr = e1,r+1 − e1,r + er+1,t (∀1 ≤ r ≤ t − 2). We t−1 t−1 claim that w1 − ⋃t−1 s=3 Ws , w2 − ⋃s=3 Ws , . . . , wt−2 − ⋃s=3 Ws are linearly independent. Otherwise, there exist θ1 , . . . , θt−2 , not all zero, such that θ1 w1 + θ2 w2 + ⋯ + θt−2 wt−2 ∈ ⋃l−1 s=3 Ws . Note that T T for any 1 ≤ r ≤ t − 2, er+1,t (θ1 w1 + ⋯ + θt−2 wt−2 ) = θr . On the other hand, er+1<l w = 0 for any w ∈ ⋃t−1 s=3 Ws , which implies θr = 0 for every 1 ≤ r ≤ t − 2, contradiction! Hence w1 , w2 , . . . , wt−2 ̃ span a space of dimension t − 2 outside ⋃t−1 s=3 Ws , which means dim (Wt ) ≥ t − 2. t−1 ̃t . In fact, for any 1 ≤ j < k < (2) We then show that wr − ⋃s=3 Ws (1 < r ≤ t−2) form a basis of W t, ejk − ejt + ekt = −(e1j − e1t + ejt ) + (e1k − e1t + ekt ) + (e1j − e1k + ejk ) ∈ −wj−1 + wk−1 + ⋃t−1 s=3 Ws . t−1 Thus Wt = span{ejk − ejt + ekt ∶ ∀1 ≤ j < k < t} = span{w1 , . . . , wt−2 } + ⋃s=3 Ws . Hence ̃t = span{w1 , . . . , wt−2 } − ⋃t−1 W s=3 Ws , which is of dimension t − 2. m ̃t ) = ∑m In summary, dim(V ) = ∑s=t dim(W l=3 (t − 2) = (m − 1)(m − 2)/2. Note that dim(U ) + m ̃f is in either U or V. ( ) dim(V ) = (m − 1) + (m − 1)(m − 2)/2 = 2 = q, so any eigenvector of M ̃f ), the result follows. And as Mf = diag(1, M Proof of Proposition 3. For any two components j1 ≠ j2 , we write “j1 precedes j2 ” as “j1 ≺≺ j2 ”. Then for any 1 ≤ j < k < l ≤ m the transitive property implies: ND (j ≺≺ k, k ≺≺ l, j ≺≺ l) = ND (j ≺≺ k, k ≺≺ l), and ND (j ≺≺ l, l ≺≺ k, j ≺≺ k) = ND (j ≺≺ l, l ≺≺ k), 14 where for any order relationship R, ND (R) represents the number of treatments in D for which R holds. Then: ND (j ≺≺ k, j ≺≺ l) = ND (j ≺≺ k, j ≺≺ l, k ≺≺ l) + ND (j ≺≺ k, j ≺≺ l, l ≺≺ k) = ND (j ≺≺ k, j ≺≺ l, k ≺≺ l) + ND (j ≺≺ l, l ≺≺ k, j ≺≺ k) = ND (j ≺≺ k, k ≺≺ l) + ND (j ≺≺ l, l ≺≺ k). Similarly, ND (k ≺≺ j, l ≺≺ j) = ND (k ≺≺ j, l ≺≺ k) + ND (l ≺≺ j, k ≺≺ l), and therefore ND (j ≺≺ k, j ≺≺ l) + ND (k ≺≺ j, l ≺≺ j) = [ND (j ≺≺ k, k ≺≺ l) + ND (k ≺≺ j, l ≺≺ k)] + [ND (j ≺≺ l, l ≺≺ k) + ND (l ≺≺ j, k ≺≺ l)]. (17) ̃ 1 k1 , j2 k2 ) = (2/N )[ND (j1 ≺≺ k1 , j2 ≺≺ k2 ) + ND (k1 ≺≺ Note that for any j1 < k1 , j2 < k2 , M(j j1 , k2 ≺≺ j2 )] − 1 = 1 − (2/N )[ND (j1 ≺≺ k1 , k2 ≺≺ j2 ) + ND (k1 ≺≺ j1 , j2 ≺≺ k2 )]. Hence from (17) we ̃ ̃ ̃ get M(jk, jl) − M(jk, kl) + M(jl, kl) = 1, and the result follows. T ̃ ̃ ̃ ̃ ̃ Proof of Corollary 1. For any jkl ∈ T , vjkl Mvjkl = M(jk, jk)+M(jl, jl)+M(kl, kl)−2M(jk, jl)+ ̃ ̃ ̃ ̃ ̃ 2M(jk, kl) − 2M(jl, kl). By Proposition 3, M(jk, jl) − M(jk, kl) + M(jl, kl) = 1 holds for any T ̃ PWO design. Thus vjkl Mvjkl = 1 + 1 + 1 − 2 = 1. To (computationally) prove Theorem 1 via Corollary 1, the following lemma helps. T . We have: (i) tr(V) = 3(m2 ), (ii) The eigenvalues of V are Lemma 2. Let V = ∑jkl∈T vjkl vjkl either 0 or m. T T vjkl ) = 3(m2 ). Now we study ) = ∑jkl∈T (vjkl Proof of Lemma 2. First, tr(V) = ∑jkl∈T tr(vjkl vjkl the eigenstructure of V. For any αβ ∈ S, T eαβ ) Veαβ = ∑ vjkl (vjkl jkl∈T = ∑ vαβl − ∑ vαkβ + ∑ vjαβ l>β j<α α<k<β = ∑(eαβ − eαl + eβl ) − ∑ (eαk − eαβ + ekβ ) + ∑ (ejα − ejβ + eαβ ) l>β j<α α<k<β (18) = (m − 2)eαβ + ∑ ejα − ∑ eαj − ∑ ejβ + ∑ eβj j<α j>α,j≠β j<β,j≠α = meαβ + ∑ ejα − ∑ eαj − ∑ ejβ + ∑ eβj . j<α j>α j<β j>β j>β From (18) it is evident that V(eαβ − eαγ + eβγ ) = m(eαβ − eαγ + eβγ ) for any αβγ ∈ T . Hence any vαβγ is an eigenvector of V with eigenvalue m. Again denote V = span{vαβγ ∶ αβγ ∈ T }, ). Noting that m ⋅ (m−1 ) = 3(m2 ) = tr(V) and that V is non-negative and recall that dim(V ) = (m−1 2 2 definite, the result follows. 15 1 wT Proof of Theorem 1 via Corollary 1. Suppose M = [ ̃ ]. Instead of all the concave and w M signed permutation invariant criteria, here we only consider the D-, A-, E-, M.S.-criteria. We show that the Mf is the optimal under the four criteria, and that the D-, A-, and M.S.-optimalities T ̃ T −1 T ̃ T −1 T ∣, M−1 = [1 + w (M − ww ) w −w (M − ww ) ], ̃ hold if and only if M = Mf . Note that ∣M∣ = ∣M−ww ̃ − wwT )−1 w ̃ − wwT )−1 −(M (M ̃ 1 + wT w wT + wT M −1 2 ̃ ̃−1 and M2 = [ ̃ ̃2 ], we have ∣M∣ ≥ ∣M∣, tr (M ) ≤ 1 + tr(M ), tr (M ) ≥ w + Mw wwT + M ̃2 ), with the equalities holding if and only if w = 0. Now it suffices to study when M ̃ 1 + tr(M attains the optimum. ̃ is M ̃ = ∑q λi xi xT , where λ1 ≤ λ2 ≤ ⋯ ≤ λq , and xi ’s Suppose the eigen-decomposition of M i=1 i q T are ortho-normal. By Corollary 1, 1 = vjkl ∑i=1 λi xi xTi vjkl for every jkl ∈ T , which implies that T (m3 ) = ∑qi=1 λi (xTi ∑jkl∈T (vjkl vjkl )xi ). Denote θi = xTi Vxi . Lemma 2 indicates that (i) 0 ≤ θi ≤ m for any i, and (ii) ∑qi=1 θi = ∑qi=1 tr(xTi Vxi ) = tr(V ∑qi=1 xi xTi ) = tr(V) = 3(m3 ). In summary, we have: q m ∑ λi θi = ( ), 3 i=1 q m ∑ θi = 3( ), and 3 i=1 0 ≤ θi ≤ m for any 1 ≤ i ≤ q. (19) Denote q1 ∶= (m − 1)(m − 2)/2. Note that q1 q q1 q i=1 i=1 i=1 q1 i=q1 +1 q i=1 i=q1 +1 m ∑ λi ⋅ m − ∑ λi θi = ∑(m − θi )λi − ∑ θi λi ≤ ∑(m − θi )λq1 − ∑ θi λq1 = λq1 × (m ⋅ q1 − ∑ θi ) i=1 = λq1 ⋅ [ m m(m − 1)(m − 2) − 3( )] = 0. 3 2 1 Therefore ∑qi=1 λi ≤ (1/m) ∑qi=1 λi θi = (m3 )/m = q1 /3, and the equality holds if and only if θi = 0 ̃ for any i > q1 . Note that ∑m i=1 λi = tr(M), which obviously equals q. In summary, we have: 1 q1 1 ∑ λi ≤ , and q1 i=1 3 (20) 1 q ∑ λi = 1. q i=1 (21) 16 ̃ is no greater than 1 . Now Inequality (20) immediately implies that the minimum eigenvalue of M 3 ̃2 ) = ∑q λ2 . First, ̃−1 ) = ∑q λ−1 , and (iii) tr(M ̃ = ∏q λi , (ii) tr(M we study (i) ∣M∣ i=1 i i=1 i i=1 1 q ∑i=q1 +1 λi ∑ 1 λi ) ∏ λi = ∏ λi ∏ λi ≤ ( i=1 ) ( q1 q − q1 i=1 i=1 i=q1 +1 q q1 q q q−q1 q (22) 1 Let x = ∑qi=1 λi , then by (21) and (22), ∏qi=1 λi ≤ f (x) = (x/q1 ) 1 [(q − x)/(q − q1 )] 1 . Taking derivative one can see that f (x) is an increasing function as x ∈ (0, q1 ). Then by (20), f (x) ≤ f (q1 /3). Hence q q−q 1 1 ̃ ≤ ( q1 /3 ) ( q − q1 /3 ) ∣M∣ q1 q − q1 1 m + 1 m−1 (m + 1)m−1 ) = . = q1 ( 3 3 3q q q−q The equality holds if and only if (a) θi = 0 for any i > q1 , (b) λ1 = λ2 = ⋯ = λq1 , (c) λq1 +1 = q 2 λq1 +2 = ⋯ = λq . Similarly, we can show that ∑qi=1 λ−1 i and ∑i=1 λi attain the optimum if and only if aforementioned conditions (a), (b), and (c) are satisfied. ̃=M ̃f . The three condiIt suffices to verify that the condition (a)+(b)+(c) is equivalent to M ̃=M ̃f . On the other hand, suppose the three conditions hold. tions are obviously satisfied if M 2 T xi ) for any i > q1 . Condition (a) implies that 0 = θi = xTi Vxi = ∑jkl∈T xTi vjkl viT xi = ∑jkl∈T (vjkl Therefore for any i > q1 , xi is orthogonal to all vjkl ’s (hence the eigenspace V ). And by Proposition 2, xi (i > q1 ) must belong to the other eigenspace U . Since dim(U ) = m − 1 = q − q1 , xi ’s (i > q1 ) span the whole space U . Thus xi ’s (i ≤ q1 ) span the orthogonal complement of U , which is V . On 1 the other hand, (19) and condition (a) imply that θi = m for all i ≤ q1 . Therefore ∑qi=1 λi ⋅ m = (m3 ), 1 i.e., ∑qi=1 λi = q1 /3. Thus ∑qi=q1 +1 λi = q − q1 /3 = (m + 1)(m − 1)/3 = (m + 1)/3 ⋅ (m − q1 ). Then by condition (b) and (c), λi = 1/3 for all i ≤ q1 , and λi = (m + 1)/3 for all i > q1 . Since ̃ = 1 PV + m+1 PU . Then by V = span {xi ∶ i ≤ q1 } and U = span {xi ∶ i > q1 }, we have M 3 3 ̃=M ̃f . This completes the proof. Proposition (2), M Proof of Theorem 3. When m = 4 and r = 2, one can verify the result straightforwardly. Now consider any r ≥ 3 and m = 2r, we show that the design D given by equation 13 achieves the same information matrix as the full PWO design. First check whether each PWO factor Zij is balanced (in the PWO design of D). Without loss of generality assume jk = 12. Note that D can be divided ⎡ (m,r) (m,r) ⎤ ⎥ ⎢ Bi Bi m ⎥. For any i-th block, there are 3 cases ⎢ into L = ( r )/2 blocks, each block being ⎢ (m,r) (m,r) ⎥ ⎥ ⎢⊖Bi Bi ⎦ ⎣ (m,r) concerning the PWO factor Z12 : (i) Both component 1 and 2 appear in Bi (abbreviated as Bi ), (ii) 1 appears in Bi and 2 appears in Bi , and (iii) 2 appears in Bi and 1 appears in Bi . In every case, apparently the i-th block includes half of runs at which 1 precedes 2, i.e., Z1 2 = +. Therefore Z12 is balanced in D. Then we consider any two different PWO factors, Zj1 k1 and Zj2 k2 . First assume Zj1 k1 and Zj2 k2 are synergistic (see Section 4). Without loss of generality, assume j1 k1 = 12 and j2 k2 = 13. For the i-th block, 17 (1) If 1, 2 and 3 all appear in Bi , then apparently Z12 and Z13 are orthogonal in this block, and there are a quarter of runs (i.e., r!/2 runs) at which Z12 = Z13 = +. (2) Similarly, if 1, 2 and 3 all appear in Bi , there are r!/2 runs at which Z12 = Z13 = +. (3) If 1 appears in Bi while 2 and 3 appear in Bi . then Z12 = Z13 = + holds for all the runs (m,r) (m,r) ]), and holds for none run in the second half. Hence in the first half of block (i.e., [Bi , Bi such block contains r! runs where Z12 = Z13 = +. (4) Similarly, if 2 and 3 appear in Bi while 1 appears in Bi , the block includes r! runs at which Z12 = Z13 = +. (5) If 2 appears in Bi while 1 and 3 appear in Bi , Z12 = Z13 = + holds for none run in the first half of block, and r!/2 runs in the second half of block. (6) Similarly, Z12 = Z13 = + holds for r!/2 runs in the blocks where: (i) 1 and 3 appear in Bi while 2 appears in Bi , (ii) 3 appears in Bi while 1 and 2 appear in Bi , and (iii) 1 and 2 appear in Bi while 3 appears in Bi . Now we count the number of blocks in each of the above cases. First notice a fact that (m,r) (m,r) (m,r) (B1 , B2 , . . . , BL ) is a rearrangement of B(m,r) ’s with L + 1 ≤ i ≤ (mr) (Recall Defii nition 1, 2, 3 and that L = (mr)/2). ) Among all the Bi ’s (1 ≤ i ≤ (mr)), there are (m−3 r−3 ones that include all of component 1, 2 and 3. From the aforementioned fact, any such Bi falls in either {Bi ∶ 1 ≤ i ≤ L} or {Bi ∶ 1 ≤ i ≤ L}. As ). The aforementioned fact a result, the numbers of blocks in case (1) and case (2) add up to (m−3 r−3 m−3 also implies: the total number of blocks in case (3) and (4) is ( r−1 ); the total number of blocks in ); the total number of blocks in case (6-ii) and (6-iii) is also (m−3 ). Hence case (5) and (6-i) is (m−3 r−1 r−1 m−3 m−3 m−3 m−3 there are N1 = ( r−3 ) ⋅ r!/2 + ( r−1 ) ⋅ r! + ( r−1 ) ⋅ r!/2 + ( r−1 ) ⋅ r!/2 runs for which Z12 = Z13 = +. As m = 2r, one can verify that N1 = m!/3. Also note that every factor is balanced, it follows that the correlation between Z12 and Z13 is 1/3. Same argument, the correlation is −1/3 between any pair of antagonistic PWO factors. Last we show that any two PWO factors are uncorrelated if they do not share a common component; without loss of generality, consider Z12 and Z34 . For any i-th block: (1) If all of component 1, 2, 3 and 4 appear in Bi , or in Bi , then apparently Z12 and Z34 are orthogonal within this block. (2) If 1, 2, 3 appear in Bi while 4 appears in Bi , then Z12 = Z34 = + holds for r!/2 runs in [Bi , Bi ], and for none run in [ ⊖ Bi , Bi ]. Thus Z12 and Z34 are orthogonal within this block. (3) Same argument, Z12 and Z34 are orthogonal within any block for which exactly one component (among 1, 2, 3, 4) appears in Bi , or exactly one component appears in Bi . (4) Suppose 1 and 2 appear in Bi while 3 and 4 appear in Bi . Then for any 1 ≤ l ≤ r!, if Z12 ⋅ Z34 = +/− at the l-th run of [Bi , Bi ], then Z12 ⋅ Z34 = −/+ at the l-th run of [ ⊖ Bi , Bi ] (since ⊖Bi reverses every row of Bi ). As such, the inner product between Z12 and Z34 equals 0 within the whole block. Similarly, Z12 and Z34 are uncorrelated in the blocks that 3 and 4 appear in Bi while 1 and 2 appear in Bi . (5) It is easy to see that the correlation between Z12 and Z34 equals 1, if (i) 1 and 3 appear in Bi while 2 and 4 appear in Bi , or (ii) 2 and 4 appear in Bi while 1 and 3 appear in Bi . (6) It is easy to see that the correlation between Z12 and Z34 equals 1, if (i) 1 and 4 appear in 18 Bi while 2 and 3 appear in Bi , or (ii) 2 and 3 appear in Bi while 1 and 4 appear in Bi . Apparently the number of blocks in case (5) is the same as that in case (6). Thus the correlations in case (5) and (6) cancel out. Thus overall, Z12 and Z34 are uncorrelated in D. 19