aa r X i v : . [ s t a t . M E ] J un Generation of Fractional Factorial Designs
Roberto Fontana a Giovanni Pistone a a DIMAT Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy
Abstract
The joint use of counting functions, Hilbert basis and Markov basis allows to definea procedure to generate all the fractions that satisfy a given set of constraints interms of orthogonality. The general case of mixed level designs, without restrictionson the number of levels of each factor (like primes or power of primes) is studied. Thisnew methodology has been experimented on some significant classes of fractionalfactorial designs, including mixed level orthogonal arrays.
Key words:
Design of Experiments, Hilbert basis, Markov basis, Algebraicstatistics, Indicator polynomial, Counting function.
All the fractional factorial designs that satisfy a set of conditions in termsof orthogonality between factors have been described as the zero-set of asystem of polynomial equations in which the indeterminates are the com-plex coefficients of their counting polynomial functions (Pistone and Rogantin(2008), Fontana et al. (2000)). A short review of this theory can be found inFontana and Rogantin (2008). In Section 2 we report a part of it to facilitatethe reader. In Section 3 we write the problem of finding fractional factorial de-signs that satisfy a set of conditions as a system of linear equations in whichthe indeterminates are positive integers. In section 4, using 4ti2 (4ti2 team(2007)) we find all the generators of some classes of fractional factorial de-signs, including mixed level orthogonal arrays and sudoku designs. Finally, insection 5 we consider the moves between different fractions as integer valuedfunctions defined over the full factorial design. We build a procedure to movebetween fractions that use Markov basis.
Notation and background
We adopt the notation used in Pistone and Rogantin (2008) and denote: • by D j a factor with n j levels coded with the n j -th roots of the unity: D j = { ω , . . . , ω n j − } ω k = exp i πn j k ! ; • by D the full factorial design with complex coding D = D × · · · D j · · · × D m . • by D the cardinality of D . • by L the full factorial design with integer coding L = Z n × · · · × Z n j · · · × Z n m , • by α an element of Lα = ( α , . . . , α m ) α j = 0 , . . . , n j − , j = 1 , . . . , m . • by [ α − β ] the m -tuple made by the componentwise difference (cid:16) [ α − β ] n , . . . , [ α j − β j ] n j , . . . , [ α m − β m ] n m (cid:17) ;the computation of the j -th element is in the ring Z n j . • by X j the j -th component function, which maps a point to its i -th compo-nent: X j : D ∋ ( ζ , . . . , ζ m ) ζ j ∈ D j ;the function X j is called simple term or, by abuse of terminology, factor . • by X α the interaction term X α · · · X α m m , i.e. the function X α : D ∋ ( ζ , . . . , ζ m ) ζ α · · · ζ α m m ;We notice that L is both the full factorial design with integer coding and theexponent set of all the simple factors and interaction terms and α is botha treatment combination in the integer coding and a multi-exponent of aninteraction term.The full factorial design in complex coding is identified as the zero-set in C m of the system of polynomial equations X n j j − , j = 1 , . . . , m . (1)2 efinition 2.1 (1) A response f on the design D is a C -valued polynomialfunction defined on D .(2) The mean value on D of a response f , denoted by E D ( f ) , is: E D ( f ) = 1 D X ζ ∈D f ( ζ ) . (3) A response f is centered on D if E D ( f ) = 0 . Two responses f and g are orthogonal on D if E D ( f g ) = 0 , where g is the complex conjugate of g . It should be noticed that the set of all the responses is a complex Hilbert spacewith the Hermitian product: f · g = E D ( f g ) . Moreover(1) X α X β = X [ α − β ] ;(2) E D ( X ) = 1, and E D ( X α ) = 0 for α = 0.The set of functions { X α , α ∈ L } is an orthonormal basis of the complexresponses on design D . In fact L = D and, from properties (i) and (ii)above, it follows that: E D ( X α X β ) = E D ( X [ α − β ] ) = α = β α = β In particular, each response f can be represented as a unique C -linear combina-tion of constant, simple and interaction terms. This representation is obtainedby repeated applications of the re-writing rules derived from Equations (1).Such a polynomial is called the normal form of f on D . In this paper we intendthat all the computation are made using the normal form. Example 2.1
Consider the full factorial design. All the monomial re-sponses on D are , X , X , X , X X , X X , X X , X X X or, equivalently, X (0 , , , X (1 , , , X (0 , , , X (0 , , , X (1 , , , X (1 , , , X (0 , , , X (1 , , and L is L = { (0 , , , (1 , , , (0 , , , (0 , , , (1 , , , (1 , , , (0 , , , (1 , , } . .2 Fractions of a full factorial design A fraction F is a multiset ( F ∗ , f ∗ ) whose underlying set of elements F ∗ iscontained in D and f ∗ is the multiplicity function f ∗ : F ∗ → N that for eachelement in F ∗ gives the number of times it belongs to the multiset F .All fractions can be obtained by adding polynomial equations, called gener-ating equations to the design equations 1, in order to restrict the number ofsolutions. Definition 2.2 If f is a response on D then its mean value on F , denotedby E F ( f ) , is E F ( f ) = 1 F X ζ ∈F f ( ζ ) where F is the total number of treatment combinations of the fraction.A response f is centered if E F ( f ) = 0 . Two responses f and g are orthogonalon F if E F ( f g ) = 0 . With the complex coding the vector orthogonality of two interaction terms X α and X β as defined before (with respect to a given Hermitian product)corresponds to the combinatorial orthogonality (all the level combinationsappear equally often in X α X β ).We consider the general case in which fractions can contain points that arereplicated. Definition 2.3
The counting function R of a fraction F is a response definedon D so that for each ζ ∈ D , R ( ζ ) equals the number of appearances of ζ in thefraction. A − valued counting function is called indicator function of a singlereplicate fraction F . We denote by c α the coefficients of the representation of R on D using the monomial basis { X α , α ∈ L } : R ( ζ ) = X α ∈ L c α X α ( ζ ) ζ ∈ D c α ∈ C . As the counting function is real valued, we have c α = c [ − α ] . We will write c in place of c ,..., . Remark 2.1
The counting function R coincides with multiplicity function f ∗ . Proposition 2.1
Let F be a fraction of a full factorial design D and R = P α ∈ L c α X α be its counting function.
1) The coefficients c α are: c α = 1 D X ζ ∈F X α ( ζ ) ; in particular, c is the ratio between the number of points of the fractionand that of the design.(2) In a fraction without replications, the coefficients c α are related accordingto: c α = X β ∈ L c β c [ α − β ] . (3) The term X α is centered on F , i.e. E F ( X α ) , if, and only if, c α = c [ − α ] = 0 . (4) The terms X α and X β are orthogonal on F , i.e. E F ( X α X β ) = 0 , if, andonly if, c [ α − β ] = 0 . Example 2.2
We consider the fraction F = { ( − , − , , ( − , , − } of the full factorial design of Example 2.1. All the monomial responses on F andtheir values on the points are ζ X X X X X X X X X X X X ( − , − ,
1) 1 − − − − − , , −
1) 1 − − − − Using Item 1 of Proposition 2.1, it is easy to compute the coefficients c α : c (0 , , = c (0 , , = c (1 , , = c (1 , , = 0 ; c (0 , , = c (1 , , = and c (1 , , = c (0 , , = − . Hence, the indicator function is F = 12 (1 − X − X X + X X X ) . From the null coefficients we see that X and X are centered and that X isorthogonal to both X and X . (cid:3) Definition 2.4
A fraction F factorially projects onto the I -factors, I ⊂{ , . . . , m } , if the projection is a multiple full factorial design, i.e. a full fac-torial design where each point appears equally often. A fraction F is a mixedorthogonal array of strength t if it factorially projects onto any I -factors with I = t . t means that, for any choice of t columns of the matrix design, allpossible combinations of symbols appear equally often. Proposition 2.2 (Projectivity) (1) A fraction factorially projects ontothe I -factors if, and only if, all the coefficients of the counting functioninvolving only the I -factors are 0.(2) If there exists a subset J of { , . . . , m } such that the J -factors appear inall the non null elements of the counting function, the fraction factoriallyprojects onto the I -factors , with I = J c .(3) A fraction is an orthogonal array of strength t if, and only if, all thecoefficients of the counting function up to the order t are zero: c α = 0 for all α of order up to t, α = (0 , , . . . , . Example 2.3 (Orthogonal array)
The fraction of a full factorial design F O = { ( − , − , − , − , − , , ( − , − , − , , , , ( − , − , , − , − , − , ( − , − , , , , − , ( − , , − , − , − , − , ( − , , − , , , − , ( − , , , − , , , ( − , , , , − , , (1 , − , − , − , , , (1 , − , − , , − , , (1 , − , , − , , − , (1 , − , , , − , − , (1 , , − , − , , − , (1 , , − , , − , − , (1 , , , − , − , , (1 , , , , , } is an orthogonal array of strength 2; in fact, its indicator function F = 14 + 14 X X X − X X X + 18 X X X X + 18 X X X X + 18 X X X X + 18 X X X X X + 18 X X X X X + 18 X X X X X − X X X X X X contains only terms of order greater than 2, together with the constant term. (cid:3) From Proposition 2.1 and Proposition 2.2 we have that the problem of find-ing fractional factorial designs that satisfy a set of conditions in terms oforthogonality between factors can be written as a polynomial system in whichthe indeterminates are the complex coefficients c α of the counting polynomialfraction. 6 xample 3.1 Let’s consider 3 factors, each one with two levels. The indicatorfunctions F = P α c α X α such that the terms X , X , X are centered on F andthe terms X i , X j i, j = 1 , , , i = j are orthogonal on F , where F = { ζ ∈ D : F ( ζ ) = 1 } , are those for which the following conditions on the coefficients of F holds c = c + c c = 2 c c Apart from the trivial F = 0 , i.e. F = ∅ and F = 1 , i.e. F = D we find F = (1 + X X X ) and F = (1 − X X X )Let’s now introduce a different way to describe the full factorial design D andall its subsets. Let’s consider the indicator functions 1 ζ of all the single pointsof D ζ : D ∋ ( ζ , . . . , ζ m ) ζ = ( ζ , . . . , ζ m )0 ζ = ( ζ , . . . , ζ m )It follows that the counting function R of a fraction F can be written as X ζ ∈D y ζ ζ with y ζ ≡ R ( ζ ) ∈ { , , . . . , n, . . . } . The particular case in which R is anindicator function corresponds to y ζ ∈ { , } .The coefficients y ζ are related to the coefficients c α as in the following Propo-sition 3.1 Proposition 3.1
Let F be a fraction of D . Its counting fraction R can beexpressed both as R = P α c α X α and R = P ζ ∈D y ζ ζ . The relation between thecoefficients c α and y ζ is c α = 1 D X ζ ∈D y ζ X α ( ζ ) Proof.
From Proposition 2.1 we have c α = 1 D X ζ ∈F X α ( ζ ) == 1 D X ζ ∈D y ζ X α ( ζ ) ✷ .1 Strata As described in Section 2, we consider m factors, D , . . . , D m where D j ≡ Ω n j = { ω , . . . , ω n j − } , for j = 1 , . . . , m . From Pistone and Rogantin (2008),we recall two basic properties which hold true for the full design D Proposition 3.2
Let X j the simple term with level set Ω n j = { ω , . . . , ω n j − } .Let’s consider the term X rj and let’s define s j = r = 0 n j /gcd ( r, n j ) r > Over D , the term X rj takes all the values of Ω s j equally often. Proposition 3.3
Let X α = X α · · · X α m m an interaction. X α i i takes values in Ω s i where s i is determined according to the previous Proposition 3.2. Let’sdefine s = lcm ( s , . . . , s m ) . Over D , the term X α takes all the values of Ω s equally often. Let’s now define the strata that are associated to simple and interaction terms.
Definition 3.1
Given a term X α , α ∈ L = Z n × . . . × Z n m the full design D is partitioned into the the following strata D αh = n ζ ∈ D : X α ( ζ ) = ω h o where ω h ∈ Ω s and s is determined according to the previous Propositions 3.2and 3.3. Remark 3.1
We define strata using the conjugate X α of the term in place ofthe term X α itself because it will simplify the notations. Remark 3.2
Each stratum is a regular fraction whose defining equation is X α ( ζ ) = ω − h , Pistone and Rogantin (2008). We use n α,h to denote the number of points of the fraction F that are in thestratum D αh , with h = 0 , . . . , s − n α,h = X ζ ∈ D αh y ζ The following Proposition 3.4 links the coefficients c α with n α,h . Proposition 3.4
Let F be a fraction of D with counting fraction R = P α ∈ L c α X α . ach c α , α ∈ L , depends on n α,h , h = 0 , . . . , s − , as c α = 1 D s − X h =0 n α,h ω h where s is determined by X α (see Proposition 3.3). Viceversa, each n α,h , h =0 , . . . , s − , depends on c [ − kα ] , k = 0 , . . . , s − as n α,h = D s s − X k =0 c [ − kα ] ω [ hk ] Proof.
Using Proposition 3.1, it follows that we can write the coefficients c α in the following way c α = 1 D X ζ ∈D y ζ X α ( ζ ) = 1 D s − X h =0 ω h X ζ ∈ D αh y ζ = 1 D s − X h =0 n α,h ω h For the viceversa, we observe the indicator function of strata can be obtainedas follows. We define ˜ F s ( ζ ) = s − X k =0 ζ k = − ζ s − ζ if ζ = 1 s if ζ = 1We have ˜ F s ( ω k ) = 0 for all ω k ∈ Ω s , k = 0. It follows that F α, ( ζ ) = 1 s ˜ F s ( ζ α ) = 1 s (cid:16) ζ α + . . . + ζ ( s − α (cid:17) is the indicator function associated to D α .The indicator of D αh = n ζ ∈ D : X α ( ζ ) = ω h o = n ζ ∈ D : X α ( ζ ) = ω [ − h ] o willbe F α,h ( ζ ) = F s ( ω h ζ α ) = 1 s (cid:16) ω h ζ α + . . . + ω [( s − h ] ζ ( s − α (cid:17) We get n α,h = X ζ ∈ D αh R ( ζ ) = X ζ ∈D F α,h ( ζ ) R ( ζ ) == X ζ ∈ D s s − X k =0 ω [ kh ] X kα ( ζ ) ! X β c β X β ( ζ ) == D s X k,β :[ kα + β ]=0 ω [ kh ] c β = D s s − X k =0 ω [ kh ] c [ − kα ] ✷ emark 3.3 From Proposition 3.4 we get n ,h = 0 , h = 1 , . . . , s − n α, = D s s − X k =0 c [ − kα ] and in particular n , = F . We now use a part of Proposition 3 of Pistone and Rogantin (2008) to getconditions on n α,h that makes X α centered on the fraction F . Proposition 3.5
Let X α be a term with level set Ω s on full design D . Let P ( ζ ) the complex polynomial associated to the sequence ( n α,h ) h =0 ,...,s − so that P ( ζ ) = s − X h =0 n α,h ζ h and let’s denote by Φ s the cyclotomic polynomial of the s -roots of the unity.(1) Let s be prime. The term X α is centered on the fraction F if, and onlyif, its s levels appear equally often: n α, = n α, = . . . = n α,s − = λ α (2) Let s = p h . . . p h d d with p i prime, for i = 1 , . . . , d . The term X α is centeredon the fraction F if, and only if, the remainder H ( ζ ) = P ( ζ ) mod Φ s ( ζ ) whose coefficients are integer linear combinations of n α,h , h = 0 , . . . , s − ,is identically zero.Proof. See Proposition 3 of Pistone and Rogantin (2008). ✷ Remark 3.4
Being D αh a partition of D , if s is prime we get λ α = F s . If we remind that n α,h are related to the values of the counting function R ofa fraction F by the following relation n α,h = X ζ ∈ D αh y ζ , this Proposition 3.5 allows to express the condition X α is centered on F asinteger linear combinations of the values R ( ζ ) of the counting function overthe full design D . In the Section 4, we will show the use of this property togenerate fractional factorial designs. 10e conclude this section limiting to the particular case where all factors havethe same number of levels s and s is prime. We provide some results concerningthe coefficients of counting functions, regular fractions, wordlength patternsand margins. From Proposition 3.5 we get the following result on the coefficients of a count-ing function
Proposition 3.6
Given a counting function R = P α c α X α , if c α = 0 then c [ k · α ] = 0 for all k = 1 , . . . , s − , where [ k · α ] is α + . . . + α | {z } k times in the ring Z ms .Proof. Let’s consider c k · α . From Proposition 3.5, c k · α is equal to zero if, andonly if, X ζ ∈ D k · α y ζ = X ζ ∈ D k · α y ζ = . . . = X ζ ∈ D k · αs − y ζ We observe that D k · αh = n ζ ∈ D : X k · α ( ζ ) = ω h o == (cid:26) ζ ∈ D : X α ( ζ ) k = ω h (cid:27) = n ζ ∈ D : X α ( ζ ) = ω [ kh ] o = D α [ kh ] where [ kh ] is h + . . . + h | {z } k times in the ring Z s .It follows that X α and X k · α partition D in the same strata and therefore weget the proof. ✷ Let’s consider a fraction F without replicates and with indicator function F = P α c α X α . Proposition 5 in (Pistone and Rogantin (2008)) states that afraction F is regular if, and only if, its indicator function F has the form F = 1 l X α ∈L e ( α ) X α where L ⊆ L , L is a subgroup of L and e : L → { ω , . . . , ω s − } is a givenmapping. 11f we use Proposition 3.5 we immediately get a characterisation of regularfractions based on the frequencies n α,h . Proposition 3.7
Given a single replicate fraction F with indicator function F = P α c α X α the following statements are equivalent:(i) F is regular(ii) for n α,h there are only two possibilities(a) if c α = 0 then n α,h = F s , h = 0 , . . . , s − ,(b) if c α = 0 then ∃ h ∗ ∈ { , . . . , s − } such that n α,h = D l if h = h ∗ otherwiseProof. Using Proposition 3.4 we get c α = 1 D s − X h =0 n α,h ω h Proposition 5 in Pistone and Rogantin (2008) gives the following conditionson the coefficients of the indicator function F of a regular fraction F : c α = e ( α ) l , α ∈ L ⊆ L e : L → { ω , . . . , ω s − } , l = mathcalL and L is a subgroup of L .Let’s consider α ∈ L . We get1 D s − X h =0 n α,h ω h = e ( α ) l Let’s suppose e ( α ) = ω h ∗ . We obtain1 D s − X h =0 ,h = h ∗ n α,h ω h + ( 1 D n α,h ∗ − l ) ω h ∗ = 0 (2)To simplify the notation we let a h = D n α,h , h = 0 , . . . , s − , h = h ∗ and a h ∗ = D n α,h ∗ − l . Therefore, from the proof of item (1) of Proposition 3.5,for the relation 2 to be valid, it should be a = a = . . . = a s − Being P s − h =0 n α,h = F it follows s − X h =0 n α,h = s − X h =0 ,h = h ∗ ( D ) a h + ( D )( a h ∗ + 1 l ) = ( D ) s − X h =0 a h + ( D ) l = F a h = 1 s ( D ) ( F − ( D ) l )We finally get n α,h = s ( F − ( D ) l ) + ( D ) l if h = h ∗ s ( F − ( D ) l ) otherwiseBeing L a subgroup of L it follows that 0 ∈ L and so c = 1 /l . We also knowthat c = F D and therefore F = D l For the null coefficients of F , { c α : α ∈ L − L} , it is enough to use Proposition ?? to conclude the proof. ✷ Aberration is often used as a criterion to compare fractional factorial de-signs. The generalized minimum aberration, proposed by Xu and Wu (2001),is based on the generalised wordlength pattern, see also Beder and Willenbring(2009). It can be shown that the generalized wordlengths can be written interms of the squares of the modules of the coefficients c α , obtaining A j = D F ! X wt ( α )= j | c α | = 1 c X wt ( α )= j | c α | for j = 1 , . . . , m where wt ( α ) is the Hamming weight of α , i.e. the number of nonzero compo-nents of α . We now express the square of the module of the coefficient c α interms of n α,h . Proposition 3.8 | c α | = 1( D ) s − X h =0 ( n α,h − n α,h n [ α,h − γ ] ) for γ ∈ { , . . . , s − } Proof.
From Proposition 3.4 we get c α = 1 D s − X h =0 n α,h ω h It follows 13 c α | = c α c α == 1( D ) ( s − X h =0 n α,h ω h )( s − X k =0 n α,k ω k ) == 1( D ) ( s − X h =0 n α,h ω h )( s − X k =0 n α,k ω [ s − k ] ) == 1( D ) s − X γ =0 s − X p =0 n α,p n [ α,p − γ ] ω γ | c α | must be a real number. Being ω = 1 it follows( 1( D ) s − X p =0 n α,p − | c α | ) ω + 1( D ) s − X γ =1 s − X p =0 n α,p n [ α,p − γ ] ω γ = 0 (3)To simplify the notation we let a = ( D ) P s − p =0 n α,p − | c α | ) and a γ = D ) P s − p =0 n α,p n [ α,p − γ ] , γ = 1 , . . . , s −
1. Therefore, by Lemma ?? , for the re-lation 3 to be valid, it should be a = a = . . . = a s − Using one of the equalities, a = a h h = 1 , . . . , s −
1, it follows | c α | = 1( D ) s − X p =0 ( n α,p − n α,p n [ α,p − h ] ) ✷ Remark 3.5
Proposition 3.8 provides a useful tool to compute the modulesof the coefficients c α . Indeed it is enough to choose γ = 1 and compute | c α | as D ) P s − h =0 ( n α,h − n α,h n [ α,h − ) ; Remark 3.6
We make explicit these relations for and level fraction.If s = 2 then | c α | = 1( D ) ( n α, − n α, ) If s = 3 then, choosing γ = 1 , | c α | = 1( D ) ( n α, + n α, + n α, − n α, n α, − n α, n α, − n α, n α, ) Remark 3.7
We observe that, denoting by n α the mean of the values of n α,h , α = s P s − h =0 n α,h , we get s − X h =0 ( n α,h − n α ) = s − X h =0 n α,h − sn α We have n α = 1 s s − X h,k =0 n α,h n α,k == 1 s s − X h =0 n α,h + 2 s − X h =0 n α,h n α, [ h − + . . . s − X h =0 n α,h n α, [ h − s ∗ ] ! where s ∗ = s − . Proposition 3.8 states that all the quantities P s − h =0 n α,h n α, [ h − γ ] are equal and so, choosing, without loss of generality, γ = 1 , we get n α = 1 s s − X h =0 n α,h + 2 s ∗ s − X h =0 n α,h n α, [ h − ! = 1 s s − X h =0 n α,h + ( s − s − X h =0 n α,h n α, [ h − ! and therefore s − X h =0 ( n α,h − n α ) = s − X h =0 n α,h − sn α == s − s s − X h =0 n α,h − s − X h =0 n α,h n α, [ h − ! == s − s ( D ) | c α | It follows that, if we denote by σ α the variance of n α,h , σ α = s P s − h =0 ( n α,h − n α ) we get | c α | = s ( s − D ) ! σ α and so the square of the module of c α represents, apart from a multiplicativeconstant, the variance of the frequencies n α,h .3.5 Margins We now examine the relationship between the margins and the coefficientsof the counting functions. We refer to (Pistone and Rogantin (2008)) and wereport here a part of it. 15or each point ζ ∈ D we consider the decomposition ζ = ( ζ I , ζ J ) where I ⊆{ , . . . , m } and J = { , . . . , m } − I ≡ I c is its complement. We denote by R I ( ζ I ) the number of points in F whose projection on the I factors is ζ I .In particular if I = { , . . . , m } we have R I = R and if I = ∅ we have R I = F .We denote by L I the subset of the exponents restricted to the I factors andby α I an element of L I : L I = { a I = ( α , . . . , α m ) , α j = 0 if j ∈ J } Then for each α ∈ L and ζ ∈ D we have α = α I + α J and X α ( ζ ) = X αI ( ζ I ) X αj ( ζ J ). Finally we denote by D I and D J the full factorial over the I factors and J factors, respectively ( D = D I × D J ).We have the following proposition (see item 1 and 2 of Proposition 4 ofPistone and Rogantin (2008)) Proposition 3.9
Given a fraction F of D (1) the number of replicates of the points of F projected on the I factors is: R I ( ζ I ) = D J X α I c α I X α I ( ζ I ) (2) F fully projects on the I factors if, and only if, R I ( ζ I ) = D J · c = D J F D = F D I We will refer to R I as k -margin, where k = I . The number of k -margins is (cid:16) mk (cid:17) and each k -margin can be computed over s k points ζ I ∈ D I . It followsthat there are (1 + s ) m marginal values in total.Using item 1 of Proposition 3.9 and reminding that we work with a primenumber of level s we have R I ( ζ I ) = s m − k X α I c α I ζ α I I or, by the definition of R I as the restriction of R over the I factors, X ζ J ∈D J R ( ζ I , ζ J ) ≡ X ζ J ∈D J y ζ I ,ζ J = s m − k X α I c α I ζ α I I We point out the following relationship between margins.
Proposition 3.10 If A ⊆ B ⊆ { , . . . , m } and R B ( ζ B ) = s m − k B c then R A ( ζ A ) = s m − k A c where B = k B and A = k A roof. Let’s put A = B − A . We have R A ( ζ A ) = X ζ A ∈ A R A ∪ A ( ζ A , ζ A ) = X ζ A ∈ A R B ( ζ A , ζ A ) = s k B − k A s m − k B c = s m − k A c ✷ We finally observe that, as we already pointed out, given
C ⊆ L a set ofconditions c α = 0 , α ∈ C translates in a set of conditions P ζ ∈ D αh y ζ = λ, h =0 , . . . , s − , α ∈ C where λ does not depend by α (and by h ). In general,with respect to margins, the situation is different. For example let’s supposeto have a F that fully projects over the I and the I factors, with I ∩ I = ∅ and I = I . From Proposition 3.9 we obtain R I ( ζ I ) = D s I and R I ( ζ I ) = D s I Let use strata to generate fractions that satisfy a given set of constrains onthe coefficients of their counting functions. Formally we give the followingdefinition
Definition 4.1
A counting function R = P α c α X α associated to F is a C -compatible counting function if its coefficients satisfy to c α = 0 , α ∈ C , C ⊆ Z n × . . . Z n m We will denote by OF ( n . . . n m , C ) the set of all the fractions whose countingfunctions are C -compatible.In the next sections, we will show our methodology on Orthogonal Arrays andSudoku designs. OA ( n, s m , t )Let’s consider OA ( n, s m , t ), i.e. orthogonal arrays with n rows and m columnswhere each columns has s symbols, s prime and with strength t .Using Proposition 2.2 we have that the coefficients of the corresponding count-ing functions must satisfy the conditions c α = 0 for all α ∈ C where C ⊆ L = { α : 0 < k α k ≤ t } where k α k is the number of non null elements of α . Wehave N = P tk =1 (cid:16) mk (cid:17) ( s − k coefficients that must be null.17t follows that OF ( s m , C ) = S n OA ( n, s m , t ).Now using Proposition 3.5, we can express these conditions using strata. If weconsider α ∈ C we write the condition c α = 0 as P ζ ∈ D α y ζ = λ P ζ ∈ D α y ζ = λ. . . P ζ ∈ D αs − y ζ = λ To obtain all the conditions it is enough to vary α ∈ C . We use Proposition 3.6to limit to the α that give different strata. It is easy to show that we obtain N = N s − different α , each of them generate s linear equations, for a total of N = sN = s t X k =1 mk ! ( s − k − constraints on the values of the counting function over D .We therefore get the following system of linear equations AY = λ A is the ( N × s m ) matrix whose rows contains the values, over D , of theindicator function of the strata, 1 D αh , Y is the s m column vector whose entriesare the values of the counting function over D , λ will be equal to F s and1 is the s m column vector whose entries are all equal to 1. We can write anequivalent homogeneous system if we consider λ as a new variable. We obtain˜ A ˜ Y = 0where ˜ A = A − − . . . − = [ A, − Y = Yλ = ( Y, λ )In an equivalent way, we can also express the conditions c α = 0 for all α ∈ C
18n terms of margins. We obtain R I ( ζ I ) = s m − ( I ) c where I ⊆ { , . . . , m } and 1 ≤ I ≤ t . If we recall Proposition 3.10, we canlimit to the margins R I where I = t . We have s t (cid:16) mt (cid:17) values of such t margin X ζ J ∈D J y ζ I ,ζ J = s m − t c In this case, with the same approach that we adopted for strata, we obtain asystem of linear equations BY = ρ ρ = s m − t c and its equivalent homogeneous system˜ B ˜ Y = 0Now we can find all the generators of OF ( s m , C , that means of OrthogonalArrays OA ( n, s m , t ), by computing the Hilbert Basis corresponding to ˜ A (or,equivalently, to ˜ B ). This approach is the same of Carlini and Pistone (2007)but, in that work, the following conditions were used c α = 1 D X ζ ∈F X α ( ζ ) = 1 D X ζ ∈D X α ( ζ ) y ζ = 0The advantage of using strata (or margins) is that we avoid computations withcomplex numbers ( X α ( ζ )). We explain this point in a couple of examples. Forthe computation we use 4ti2 (4ti2 team (2007)).We use both ˜ A (strata) and ˜ B (margins) because, even if they are fully equiv-alent from the point of view of the solutions that they generate, they performdifferently from the point of view of the computational speed. OA ( n, , OA ( n, ,
2) were investigated in Carlini and Pistone (2007). We build boththe matrix ˜ A and ˜ B . They have 30 rows and 40 rows, respectively and 33columns. We find the same 26 ,
142 solutions as in the cited paper. OA ( n, , A and ˜ B . They have 54 rows and 27 rows, respec-tively and 28 columns. We find 66 solutions, 12 have 9 points, all different and54 have 18 points, 17 different. 19inally we point out that 4ti2 allows to specify upper bounds for variables.For example, if we use ˜ B and we are interested in single replicate orthogonalarrays, we can set 1 as the upper bound for y ζ , ζ ∈ D . The upper bound forthe variable ρ can be set to s m − t ≡ − that corresponds to c = 1, i.e. to thefull design D . OA ( n, n . . . n m , t )Let’s now consider the general case in which we do not put restrictions on thenumber of levels. OA ( n, , . Using Proposition2.2 we have that the coefficients of the corresponding counting functions mustsatisfy the conditions c α = 0 for all α ∈ C where C ⊆ L = { α : k α k = 1 } .Let’s consider c , . From Proposition 3.2 we have that X takes the values inΩ s where s = 4. From Proposition 3.5, X will be centered on F if, and onlyif, the remainder H ( ζ ) = P ( ζ ) mod Φ ( ζ )is identically zero. We have Φ ( ζ ) = 1 + ζ (see Lang (1965)) and so we cancompute the remainder H ( ζ ) = n (1 , , − n (1 , , + ( n (1 , , − n (1 , , ) ζ The condition H ( ζ ) identically zero translates into n (1 , , − n (1 , , = 0 n (1 , , − n (1 , , = 0Let’s now consider c , . From Proposition 3.2 we have that X takes the valuesin Ω s where s = 2. From Proposition 3.5, X will be centered on F if, andonly if, the remainder H ( ζ ) = P ( ζ ) mod Φ ( ζ )is identically zero. We have Φ ( ζ ) = 1 + ζ (see Lang (1965)) and so we cancompute the remainder H ( ζ ) = n (2 , , − n (2 , , If we repeat the same procedure for all the α such that k α k = 1 and we recallthat n α,h = X ζ ∈ D αh y ζ OA ( n, ,
1) become the integer solutions of the followinginteger linear homogeneous system − − − − − − − − − − − − − − − −
11 0 − − − − − − − − − − − − − − − −
11 1 1 1 − − − − − − − −
11 1 1 1 0 0 0 0 − − − − − − − − y y y y y y y y y y y y y y y y Using 4ti2 we find 24 solutions that correspond to all the Latin HypercupeDesigns (LHD). OA ( n, , c α = 0 forall α ∈ C where C ⊆ L = { α : k α k = 1 } .Let’s consider c , . From Proposition 3.2 we have that X takes the values inΩ s where s = 6. From Proposition 3.5, X will be centered on F if, and onlyif, the remainder H ( ζ ) = P ( ζ ) mod Φ ( ζ )is identically zero. We have Φ ( ζ ) = 1 − ζ + ζ (see Lang (1965)) and so wecan compute the remainder H ( ζ ) = n (1 , , − n (1 , , − n (1 , , + n (1 , , + ( n (1 , , + n (1 , , − n (1 , , − n (1 , , ) ζ
21f we repeat the same procedure for all the α such that k α k = 1 and we recallthat n α,h = X ζ ∈ D αh y ζ orthogonal arrays OA ( n, ,
1) become the integer solutions of an integer linearhomogeneous system AR = 0 where the matrix A is built as in the previouscase of OA ( n, , As shown in Fontana and Rogantin (2008), a sudoku can be described usingits indicator function. Here we report a very short synthesis of Section 1.3 ofthat work.A p × p with p prime sudoku design can be seen as a fraction F of the fullfactorial design D : D = R × R × C × C × S × S where each factor is coded with the p -th roots of the unity. R and R , C and C , S and S , represent the rows, the columns and the symbols of the sudokugrid, respectively.The following proposition (Proposition 5 of Fontana and Rogantin (2008))holds. Proposition 4.1
Let F be the indicator function of a fraction F of a design design , F = P α ∈ L b α X α . The fraction F corresponds to a sudoku grid if andonly if the coefficients b α satisfy the following conditions:(1) b = 1 /p , i.e. the ratio between the number of points of the fractionand the number of points of the full factorial design is /p ;(2) for all i j ∈ { , , . . . , p − } :(a) b i i i i = 0 for ( i , i , i , i ) = (0 , , , ,(b) b i i i i = 0 for ( i , i , i , i ) = (0 , , , ,(c) b i i i i = 0 for ( i , i , i , i ) = (0 , , , ,(d) b i i i i = 0 for ( i , i , i , i ) = (0 , , , i.e. the fraction factorially projects onto the first four factors and ontoboth symbol factors and row/column/box factors, respectively. From this Proposition, we define C as the union of C , C , C and C , where22 = { ( i i i i
00) : ( i , i , i , i ) = (0 , , , }C = { ( i i i i ) : ( i , i , i , i ) = (0 , , , }C = { (00 i i i i ) : ( i , i , i , i ) = (0 , , , }C = { ( i i i i ) : ( i , i , i , i ) = (0 , , , } The problem of finding Sudoku becomes equivalent to find C -compatible count-ing functions, that are (i) indicator functions and (ii) that satisfy the additionalrequirement b = 1 /p . × Sudoku
We use the conditions C to build both the matrices ˜ A and ˜ B . ˜ A has 78 rows.With respect to ˜ B , that corresponds to the margins that must be constant,if we recall Proposition 3.10 we obtain 64 constraints, all corresponding to4-margins.To find all sudoku we use 4ti2, specifying the upper bounds for all the 65variables. The upper bounds for y ζ , ζ ∈ D must be equal to 1. If we use ˜ A ,the upper bound for λ must be set equal to F s ≡ = 8, while if we use ˜ b the upper bound for ρ must be set equal to s m − k b ≡ = 1.We find all the 288 different 4 × A the total time was 31 . B the total time was only 58 .
04 seconds on the samecomputer.If we admit counting functions with values in { , , } and F ≤
32 we find55 ,
992 solutions.
Sometimes, given a set of conditions C we are interested in picking up a solutionmore than in finding all the generators. The basic idea is to generate somehowa starting solution and then to randomly walk in the set of all the solutionsfor a certain number of steps, taking the arrival point as a new but still C -compatible counting function.Let’s use the previous results on strata to get a suitable set of moves . We willshow this procedure in the case in which all the factors have the same numberof levels s , S prime, but it can also be applied to the general case. In Section 4we have shown that counting functions must satisfy the following set of linear23quations AY = λ A corresponds to the set of conditions C written in terms of strata.It follows that if, given a C -compatible solution Y , such that AY = λ
1, wesearch for an additive move X such that A ( Y + X ) is still equal to λ
1, wehave to solve the following linear homogenous system AX = 0with X = ( x ζ ) , ζ ∈ D , x ζ ∈ Z and y ζ + x ζ ≥ ζ ∈ D . We observe thatthis set of conditions allows to determine new C -compatible solutions that givethe same λ . We know that λ = F s so this homogenous system determinesmoves that do not change the dimension of the solutions .Let’s now consider the extended homogeneous system, where ˜ A has alreadybeen defined in Section 4, ˜ A ˜ X = 0with ˜ X = (˜ x ζ ) , ζ ∈ D , ˜ x ζ ∈ Z and ˜ y ζ + ˜ x ζ ≥ ζ ∈ D .Given ˜ Y = ( Y, λ Y ), where Y is C -compatible counting function and λ Y = P ζ y ζ s , the solutions of ˜ A ˜ X = 0 determine all the other ˜ Y + ˜ X = ( Y + X, λ Y + X )such that ˜ A ( ˜ Y + ˜ X ) = 0. Y + X are C -compatible counting functions whosesizes, sλ Y + X , are, in general, different from that of Y . We use the theory of Markov basis (see for example Drton et al. (2009) whereit is also available a rich bibliography on this subject) to determine a set ofgenerators of the moves.We use the following procedure in order to randomly select a C -compatiblecounting function. We compute a Markov basis of ker( A ) using 4ti2 (4ti2 team(2007)). Once we have determined the Markov basis of ker( A ), we make arandom walk on the fiber of Y , where Y , as usual, contains the values ofthe counting function of an initial design F . The fiber is made by all the C -compatible counting functions that have the same size of F . The randow walkis done randomly choosing one move among the feasible ones, i.e. among themoves for which we do not get negative values for the new counting function.In the next paragraphs we consider moves for the cases that we have alreadystudied in Section 4. 24 .2 Orthogonal arrays5.2.1 OA ( n, , A , already built in Section 4.1.1 and give it as input to4ti2 to obtain the Markov Basis, that we denote by M . It contains 5 . M = ( x ζ ) ∈ M we define M + = max( x ζ ,
0) and M − = max( − x ζ , M = M + − M − .As an initial fraction F , we consider the eight-run regular fraction whoseindicator function R is R = 14 (1 + X X X )(1 + X X X )We obtain the set of feasible moves observing that a move M ∈ M , to befeasible, should be not negative when R is equal to zero that means(1 − R ) M − = 0We find 12 moves. Analogously an element M ∈ M such that(1 − R ) M + = 0gives a feasible move, − M . In this case we do not find any of such element.Therefore, given R , the set of feasible moves becomes M R that contains12 + 0 different moves.We randomly choose one move M R out of the 12 available ones and move to R = R + M R We run 1.000 simulations repeating the same loop, generating R i as R i = R i − + M R i − .We obtain all the 60 different 8-run fractions, each one with 8 different pointsas in Carlini and Pistone (2007).Using ˜ A we obtain the set ˜ M that contains 18 different moves. OA ( n, , A as built in the Section 4.1.2, we use 4ti2 to generate the Markovbasis corresponding to the homogeneous system AX = 0. We obtain M thatcontains 81 different moves. 25s an initial fraction we can consider the nine-run regular fraction F whoseindicator function R is R = 13 (1 + X X X + X X X )We run 1 .
000 simulations repeating the same loop, i.e. generating R i as R i = R i − + M R i − .We obtain all the 12 different 9-run fractions, each one with 9 different pointsas known in the literature and as found in Section 4.1.2.Using ˜ A we also obtain the set ˜ M that contains 10 different moves. × sudoku Using the matrix A built in Section 4.3.1, we run 4ti2 getting the Markovbasis M that contains 34 .
920 moves.We randomly choose an initial sudoku3 2 4 14 1 3 22 3 1 41 4 2 3The corresponding indicator function is F = 14 (1 − R C S S )(1 − R C S ) . Then we extract from M the feasible moves. We obtain a subset M F thatcontains 5 different moves. We repeat the procedure on −M and we obtainother 9 moves.We randomly choose one move M F out of the 5 + 9 available ones and moveto F = F + M F We run 1 .
000 simulations repeating the same loop F i = F i − + M F i − .We obtained all the 288 different 4 × Conclusions
We considered mixed level fractional factorial designs. Given the countingfunction R of a fraction F we translated the constraint c α = 0, where c α is ageneric coefficient of its polynomial representation R = P α c α X α , into a set oflinear constraints with integer coefficients on the values y ζ that R takes on allthe points ζ ∈ D . We obtained the set of generators of the solutions of someproblems using Hilbert Basis. We also studied the moves between fractions.We characterized these moves as the solution of a homogeneous linear system.We defined a procedure to randomly walk among the solutions that is based onthe Markov basis of this system. We showed the procedure on some examples.Computations have been made using 4ti2 (4ti2 team (2007)).Main advantages of the procedure are that we do not put restrictions on thenumber of levels of factors and that it is not necessary to use software thatdeals with complex polynomials.One limit is in the high computational effort that is required. In particularonly a small part of the Markov basis is used because of the requirement thatcounting functions can only take values greater than or equal to zero. Thepossibility to generate only the moves that are feasible could make the entireprocess more efficient and is part of current research. References