Finding D-optimal designs by randomised decomposition and switching

Finding D-optimal designs b y randomised decomp osition and switc hing Ric hard P . Bren t Australian National Univ ersity Dedicated to Kath y Horadam on the o ccasion of her sixtieth birthday Abstract A square { +1 , − 1 } -matrix of order n with maximal determinan t is called a satur ate d D-optimal design . W e consider some cases of satu- rated D-optimal designs where n > 2, n 6≡ 0 mo d 4, so the Hadamard b ound is not attainable, but b ounds due to Barba or Ehlich and W o j- tas may b e attainable. If R is a matrix with maximal (or conjectured maximal) determinan t, then G = R R T is the corresp onding Gr am ma- trix . F or the cases that w e consider, maximal or conjectured maximal Gram matrices are known. W e sho w ho w to generate many Hadamard equiv alence classes of solutions from a given Gram matrix G , using a randomised decomposition algorithm and ro w/column switc hing. In particular, we consider orders 26, 27 and 33, and obtain new satu- rated D-optimal designs (for order 26) and new conjectured saturated D-optimal designs (for orders 27 and 33). 1 In tro duction The Hadamar d maximal determinant (maxdet) problem is to find the maxi- m um determinan t D ( n ) of a square { +1 , − 1 } -matrix of giv en order n . Suc h a matrix A with maximal | det( A ) | is called a satur ate d D-optimal design of order n . W e are only concerned with the absolute v alue of the determinant, as the sign may b e changed b y a row or column interc hange. Hadamard [9] sho wed that D ( n ) ≤ n n/ 2 , and this b ound is attainable for n > 2 only if n ≡ 0 mo d 4. The “Hadamar d c onje ctur e” (due to P aley [20]) is that Hadamard’s b ound is attainable for all n ≡ 0 mo d 4. In this pap er w e are concerned with “non-Hadamard” cases n > 2, n 6≡ 0 mo d 4. F or such orders the Hadamard b ound is not attainable, but other upp er b ounds due to Barba [1], Ehlich [7, 8] and W o jtas [26] ma y be attainable. F or lo wer b ounds on D ( n ), see Brent and Osb orn [4], and the references giv en there. W e say that tw o { +1 , − 1 } matrices A and B are Hadamar d-e quivalent (abbreviated H-equiv alen t) if B can b e obtained from A b y a signed p er- m utation of ro ws and/or columns. If A is H-equiv alen t to B or to B T then w e say that A and B are extende d Hadamar d-e quivalent (abbreviated HT- equiv alent). Note that, if A is HT-equiv alent to B , then | det( A ) | = | det( B ) | . If A is H-equiv alent to A T then w e say that A is self-dual . W e say that a Hadamard equiv alence class is self-dual if the class contains a self-dual matrix (equiv alently , if the duals 1 of all matrices in the class are also in the class). If w e know (or conjecture) D ( n ), it is of interest to find all (or most) Hadamard equiv alence classes of { +1 , − 1 } matrices with determinan t ± D ( n ). In this pap er w e consider the orders 26, 27 and 33; similar metho ds can b e used for certain other orders. In § 2 w e consider the randomised decomp osition of (candidate) Gram matrices. Then, in § 3, w e sho w how one solution can often b e used to generate other, generally not Hadamard equiv alent, solutions via switc hing. The graph of Hadamard equiv alence classes induced b y switching is defined in § 3. W e conclude with some new results for the orders 26, 27 and 33 in §§ 4 – 6. Upp er b ounds A b ound which holds for all o dd orders, and whic h is known to b e sharp for an infinite sequence of orders ≡ 1 (mo d 4), is D ( n ) ≤ (2 n − 1) 1 / 2 ( n − 1) ( n − 1) / 2 , (1) due indep enden tly to Barba [1] and Ehlich [7]. W e call it the Barb a b ound. Brou wer [5] show ed that the Barba b ound (1) is sharp if n = q 2 + ( q + 1) 2 for q an o dd prime p o wer. The b ound is also sharp in some other cases, e.g. q = 2 and q = 4. It is not achiev able unless n is the sum of tw o consecutiv e squares. An upp er b ound due to Ehlich [8] applies only in the case n ≡ 3 (mo d 4). W e refer to [3, 8, 18, 19] for details of this b ound, which is rather complicated. The Ehlich b ound is not known to b e sharp for an y order n > 3. Another b ound, D ( n ) ≤ (2 n − 2)( n − 2) ( n − 2) / 2 , (2) 1 W e use “dual” and “transp ose” in terchangeably . 2 due to Ehlich [7] and W o jtas [26], applies in the case n ≡ 2 (mo d 4). It is kno wn to b e sharp in the following cases: (A) n = 2( q 2 + ( q + 1) 2 ), where q is an o dd prime pow er (see Whiteman [25]); and (B) n = 2( q 2 + q + 1), where q is an y (even or o dd) prime p o w er [22, 12]. Gram matrices If R is a given square matrix then the symmetric matrix G = R R T is called the Gr am matrix of R . W e may also consider the dual Gr am matrix H = R T R . Since det( G ) = det( R ) 2 , the b ounds mentioned ab o v e on det( R ) are equiv alen t to b ounds on det( G ). Indeed, this observ ation explains the form of the b ounds. F or example, the Barba b ound corresp onds to a matrix G = ( g i,j ) giv en by g i,j = n if i = j and g i,j = 1 if i 6 = j . It is easy to show via a w ell-known rank-1 up date form ula that det( G ) = (2 n − 1)( n − 1) n − 1 . Giv en a symmetric matrix G with suitable determinan t, w e say G is a c andidate Gr am matrix . It is the Gram matrix of a { +1 , − 1 } matrix if and only if it decomp oses into a pro duct of the form G = R R T , where R is a square { +1 , − 1 } matrix. 2 Decomp osition of candidate Gram matrices Supp ose that a (candidate) Gram matrix G of order n is known. W e wan t to find one or more { +1 , − 1 } matrices R such that G = RR T . Let the rows of R b e r T 1 , . . . , r T n . Then r T i r j = g i,j , 1 ≤ i, j ≤ n. If w e already know the first k rows, then we get k single-Gr am constrain ts in volving ro w k + 1: r T i r k +1 = g i,k +1 for 1 ≤ i ≤ k . These are linear constraints in the unknowns r k +1 . W e ma y b e able to find one or more solutions for row k + 1 satisfying the single-Gram constraints, or there ma y b e no solutions, in whic h case we hav e to bac ktrack. Our algorithm is describ ed in [3, § 4], so we omit the details here. W e merely note that it is p ossible to tak e adv an tage of v arious symmetries to reduce the size of the searc h space, and that it is possible to prune the search using gr am-p air constraints of the form G j +1 = RH j R T ( j > 0) if we kno w the dual Gram matrix H = R T R . In the cases considered below, there is (up 3 to signed p erm utations) only one candidate Gram matrix with the required determinan t, so there is no loss of generality in assuming that G = H . The search can b e regarded as searc hing a (large) tree with (at most) n lev els, where eac h level corresp onds to a row of R . A deterministic search t ypically searches the tree in depth-first fashion – at each no de, recursively searc h the subtrees defined by the c hildren of that no de. The aim is to find one or more lea ves at the n -th lev el of the tree, since these leav es corresp ond to complete solutions R . Deterministic, depth-first search may tak e a long time searching fruit- lessly for solutions in subtrees where no solutions exist. When G is decom- p osable, but difficult to decomp ose using a deterministic searc h, we may b e able to do b etter with a randomised search. In the randomised searc h, at each no de w e randomly c ho ose a small n umber of c hildren and recursively searc h the subtrees defined by these c hildren. A go o d choice of the a verage num b er of c hildren chosen p er no de, sa y µ , can b e determined exp erimen tally . T o o small a v alue makes it unlikely that a solution will be found; to o large a v alue mak es the search tak e to o long. W e found empirically that µ ≈ 1 . 3 works well in the cases considered b elo w. Th us, at eac h node tra versed in the searc h w e c ho ose one c hild (if there are an y) and, with probability ab out µ − 1 ≈ 0 . 3, also choose a second c hild (if there is one), then recursively searc h the subtrees defined by the selected children. F or example, in the case n = 27, there is a known Gram matrix G , due to T amura [23], whic h decomp oses into RR T , giving a { +1 , − 1 } matrix R of determinan t 546 × 6 11 × 2 26 . This determinant is conjectured to be maximal. A deterministic searc h fails to decomp ose T amura’s G in 24 hours (ex- ploring o ver 10 8 no des but reaching only depth 17 in the search tree). The tree size is probably greater than 4 × 10 9 . On the other hand, our randomised search routinely finds a decomp osi- tion of G in ab out 90 seconds. In this wa y w e ha ve found man y different H-classes of solutions. F urther details are given in § 5. Non uniformit y of sampling Unfortunately , the randomised searc h strategy describ ed ab o ve does not guaran tee that the set of H-classes of solutions (or the set of all { +1 , − 1 } matrices of the given order and determinant) is sampled uniformly . There are tw o reasons for lack of uniformity . First, the tree-generation algorithm in tro duces non-uniformit y b y taking adv antage of symmetries to reduce the size of the tree. Second, the num b er of matrices in a class is inv ersely 4 prop ortional to the order of the automorphism group of the class, so ev en if all the { +1 , − 1 } matrices were sampled uniformly , the H-classes would not necessarily b e sampled uniformly 2 . 3 Switc hing Switching is an op eration on square { +1 , − 1 } matrices which preserv es the absolute v alue of the determinan t but do es not generally preserve Hadamard equiv alence or extended Hadamard equiv alence. Th us, switching can be used to generate man y inequiv alen t maxdet solu- tions from one solution. This idea w as introduced by Denniston [6] and used to goo d effect b y Orrick [16]. Similar ideas hav e b een used by W anless [24] and others in the context of latin squares. W e only consider switching a closed quadruple of rows/columns. There are other p ossibilities, e.g. switc hing Hall sets [16]. Switc hing a closed quadruple of rows/columns Supp ose that a { +1 , − 1 } matrix R is H-equiv alent to a matrix ha ving a close d quadruple of rows, i.e. four ro ws of the form 3 :     + · · · + − · · · − − · · · − + · · · + + · · · + − · · · − + · · · + − · · · − + · · · + + · · · + − · · · − − · · · − + · · · + + · · · + + · · · + + · · · +     Then r ow switching flips the sign of the leftmost blo c k, giving     − · · · − − · · · − − · · · − + · · · + − · · · − − · · · − + · · · + − · · · − − · · · − + · · · + − · · · − − · · · − − · · · − + · · · + + · · · + + · · · +     This is H-equiv alent to flipping the signs of all but the leftmost blo c k, whic h has a nicer interpretation in terms of switching edges in the corre- sp onding bipartite graph [14, 15]. It is easy to see that ro w switc hing preserves the inner pro ducts of eac h pair of columns of R, so preserves the dual Gram matrix R T R , and hence 2 T o a certain exten t, these tw o sources of bias may tend to cancel. 3 W e write “+” for +1 and “ − ” for − 1. 5 preserv es | det( R ) | . How ever, it do es not generally preserve H-equiv alence or HT-equiv alence. Column switching is dual to row switching – instead of a closed quadruple of four ro ws, it requires a closed quadruple of four columns. Equiv alence classes generated b y switc hing, and their graphs Let A and B b e t wo H-equiv alence classes of matrices. W e say that A and B are switching-e quivalent (abbreviated “S-equiv alent”) if there exists A ∈ A and B ∈ B such that A can b e transformed to B b y a sequence of row and/or column switching op erations 4 . The size of an S-equiv alence class C , denoted b y ||C || S , is the num b er of H-equiv alence classes that it contains. If the H-equiv alence classes corresponding to matrices A and B are in the same S-equiv alence class, then w e write A ↔ B . Thus, this notation means that there is a sequence of row/column switc hes that transforms A to a matrix H-equiv alent to B . W e say that A is S-e quivalent to B . If A and B are tw o HT-equiv alence classes of matrices, then we say that A and B are ST-e quivalent if there exists A ∈ A and B ∈ B such that A can b e transformed to B by a sequence of row and/or column switc hing op erations 5 . The size of an ST-equiv alence class C , denoted b y ||C || S T , is the num b er of HT-equiv alence classes that it contains. W e say that tw o matric es A and B are ST-equiv alen t if the corresp onding HT-classes A 3 A and B 3 B are ST-equiv alent. Thus, t wo matrices A and B are ST-equiv alen t if a matrix H-equiv alent to B can b e obtained from A b y a sequence of ro w switc hes, column switches and/or transp ositions. The weight w ( H ) of a matrix H (or of the Hadamard class H 3 H ) is defined by w ( H ) = w ( H ) = 1 | Aut( H ) | , (3) where Aut( H ) is the automorphism group of H . The weight w ( C ) of an S-class C is defined b y w ( C ) = X H∈C w ( H ) . (4) The probability of finding a class by uniform random sampling of { +1 , − 1 } matrices is prop ortional to the w eight of the class, so the classes with smallest w eight are in some sense the hardest to find. (How ev er, as observed at the end of § 2, we do not sample uniformly .) 4 S-equiv alence is the same as Orric k’s Q-e quivalenc e [16, 17] in the cases that we con- sider, but the concepts are differen t if n ≡ 4 mod 8. 5 Th us, for all α ∈ A and β ∈ B , w e hav e either α ↔ β or α ↔ β T . 6 Asso ciated with an S-equiv alence class S = { H 1 , . . . , H s } of size s there is a graph 6 G = G ( S ) whose vertices are the H-classes H 1 , . . . , H s con tained in S , and where an edge connects t wo distinct vertices H i , H j if a matrix in H i can b e transformed to a matrix in H j b y a single row/column switching op eration. Similarly , for an ST-equiv alence class C = { H 1 , . . . , H s } of size s there is a graph G = G ( C ) whose v ertices are the HT-classes H 1 , . . . , H s con tained in C , and where an edge connects t wo distinct vertices H i , H j if a matrix in H i can b e transformed to a matrix in H j b y a single row/column switc hing op eration, p ossibly combined with transp osition. F or example, it is kno wn [10, 11] that there are 60 H-classes of Hadamard matrices of order 24. These form tw o S-classes, of size 1 and 59. Similarly , there are 36 HT-classes, giving tw o ST-classes, of size 1 and 35. In eac h case the class of size 1 con tains the Paley matrix, whic h has no closed quadruples. 4 Results for order 26 F or order 26 the maximal determinan t is D (26) = 150 × 6 11 × 2 25 , meeting the Ehlich-W o jtas b ound (2), and the corresp onding Gram matrix G is unique up to symmetric signed p ermutations. Without loss of generality we can assume that G has a diagonal blo c k form with blo cks of size 13 × 13 (see [7, 18, 26]). There are exactly three H-inequiv alent maxdet matrices comp osed of circulan t blo c ks [27, 13]. Ho w ever, there are man y solutions that are not comp osed of circulan t blo c ks. Orrick [17, Sec. 7] found 5026 HT-classes (9884 H-classes) of solutions by a com bination of hill-clim bing (lo cal optimisation) and switching. Using randomised decomp osition of the Gram matrix G follow ed b y switc hing, we hav e found 39 further H-classes (23 HT-classes). Th us, there are at least 9923 H-classes (5049 HT-classes) of saturated D-optimal designs of order 26. Since the randomised decomp osition program has rep eatedly found the same set of 9923 H-classes without finding any more, it is reason- able to conjecture that this is all. An exhaustiv e search to prov e this may b e feasible, but has not yet b een attempted. It is known [7, 17] that there are t wo typ es of maxdet matrices of order n = 26, related to the tw o wa ys that the row sums 2 n − 2 = 50 of the Gram matrix can b e written as a sum of squares: 50 = 7 2 + 1 2 = 5 2 + 5 2 . 6 W e ignore any lo ops or multiple edges, so all graphs considered here are simple. 7 ||C || S ||C || S T w ( C ) t yp e splits notes 8545 4323 229955/52 (5 , 5) no G = < R 3 > 7 4 9/2 (5 , 5) no new 4 3 3/2 (5 , 5) no 1 1 1/2 (5 , 5) no 1 1 1/6 (5 , 5) no new 1 1 1/78 (5 , 5) no 5, 5 5 11/6, 11/6 (5 , 5) y es 4, 4 4 2, 2 (5 , 5) y es 1, 1 1 1/2, 1/2 (5 , 5) y es new 1, 1 1 1/3, 1/3 (5 , 5) y es new 1, 1 1 1/6, 1/6 (5 , 5) y es 1310 686 3046/3 (7 , 1) no E 19 11 6 (7 , 1) no new 1 1 1/3 (7 , 1) no new 1 1 1/39 (7 , 1) no new 1 1 1/78 (7 , 1) no < R 2 > 3, 3 3 2/3, 2/3 (7 , 1) y es new 1, 1 1 1/78, 1/78 (7 , 1) y es < R 1 > 9923 5049 852013/156 — — totals T able 1: 25 S-classes and 18 ST-classes for order 26 They are called “t yp e (7 , 1)” and “t yp e (5 , 5)” respectively . The t yp e is preserv ed b y switching. If a maxdet matrix R of order 26 is normalised so that R R T = R T R = G , then λ ( R ) := X i   X j r i,j   determines the type of R : maxdet matrices of t yp e (7 , 1) hav e λ ( R ) = 182, and those of type (5 , 5) hav e λ ( R ) = 130. There are 5049 HT-classes (9923 H-classes) whic h lie in 18 ST-classes (25 S-classes). There is one “gian t” ST-class G with size ||G || S T = 4323, consisting of t yp e (5 , 5) matrices. There is another “large” ST-class E with ||E || S T = 686, consisting of t yp e (7 , 1) matrices. Eac h ST-class C of size s = ||C || S T corresp onds to either one S-class C 1 (of size ||C 1 || S < 2 s ) or tw o S-classes C 1 , C 2 (eac h of size ||C i || S = s ), dep ending on whether or not the ST-class con tains a self-dual matrix. In the former case w e say that the ST-class is self-dual , otherwise we sa y that the ST-class 8 Figure 1: The ST-class of size 11 for maxdet matrices of order 26 splits . F or example, the ST-class of size 11 is self-dual and corresponds to an S-class of size 19, but the ST-class of size 5 splits to giv e t w o S-classes of size 5. Details of all the known classes are given in T able 1. The third column of the table giv es the w eight(s) of the S-class(es) in that row, where the weigh t is defined by (4) ab o ve. The en tries lab elled “new” are not given in Orric k’s pap er [17]. The classes labelled < R i > (1 ≤ i ≤ 3), are generated b y matrices composed of circulan t blocks, using the notation of [18]. The graph asso ciated with the ST class of size 11 is shown in Figure 1. The largest automorphism group order is 22464 = 2 6 · 3 3 · 13, and all group orders divide 22464. The distribution of group orders is giv en in T able 2. In the table, the columns headed “#” give the num b er of times that the corresp onding group order o ccurs. A list of (representativ es of ) H-classes and their group orders is a v ailable from [2]. T o summarise the main results, w e hav e: Theorem 1. F or or der 26 ther e ar e at le ast 9923 Hadamar d classes of { +1 , − 1 } matric es with determinant 150 × 6 11 × 2 25 . They lie in at le ast 25 switching classes, as given in T able 1 . Pr o of. The pro of is computational. On our website [2] w e give represen- tativ es of eac h of the 18 ST-classes. F rom these “generators”, a program that implements switching can find all 5049 HT-classes; this requires only 12 iterations of ro w/column switching and taking duals. By taking duals of the 7 generators that are not self-dual, we obtain 25 generators for the 9923 H-classes. 9 Order # Order # Order # Order # 1 2823 2 4086 3 41 4 1840 6 151 8 607 12 106 16 143 24 20 32 44 36 6 39 1 48 13 64 13 72 8 78 6 96 3 108 1 156 2 216 2 288 4 576 2 22464 1 T able 2: Group orders of 9923 H-classes for order 26 5 Results for order 27 It is kno wn that the maximal determinan t D (27) for order 27 satisfies 546 ≤ D (27) 6 11 × 2 26 < 565 , where the lo wer b ound is due to T amura [23], and the upp er b ound is the (rounded up) Ehlic h b ound [8]. It is plausible to conjecture that the lo wer b ound 546 × 6 11 × 2 26 is maximal, since it is 0 . 9673 of the Ehlic h b ound and has not b een improv ed despite attempts using optimisation tec hniques that ha ve b een successful for other orders [18]. Unfortunately , proving that the lo wer b ound is maximal seems difficult – the tec hnique used in [3] to prov e analogous results for orders 19 and 37 is impractical for order 27 due to the size of the search space. T amura found a { +1 , − 1 } matrix R , with determinant 546 × 6 11 × 2 26 . The corresp onding Gram matrix G = R R T has a blo ck form with diago- nal blocks of sizes (7 , 7 , 7 , 6). Orric k [17] sho wed that T amura’s matrix R generates an ST-class T with ||T || S T = 33. The ST-class T splits in to t wo S-classes, each containing 33 H-classes. Using randomised decomp osition of T amura’s (conjectured maximal) Gram matrix, follo wed by switc hing, w e ha ve the following result. Theorem 2. Ther e ar e at le ast 6489 HT-classes (12911 H-classes ) of {± 1 } matric es of or der 27 with determinant 546 × 6 11 × 2 26 . They lie in at le ast 204 ST-classes (388 S-classes ) . The lar gest ST-class c ontains at le ast 5765 HT-classes (11483 H-classes ) . Pr o of. The pro of is computational. On our website [2] w e give represen- tativ es of eac h of the 204 ST-classes. F rom these “generators”, a program that implements switching can find all 6489 HT-classes; this requires only 28 iterations of ro w/column switching and taking duals. By taking duals of 10 ||C || S T # #split ||C || S T # #split ||C || S T # #split 5765 1 0 36 1 1 33 1 1 28 1 1 21 1 1 18 2 2 14 2 2 12 4 4 11 1 1 9 2 2 8 3 3 7 7 5 6 12 12 5 11 9 4 12 11 3 18 17 2 38 33 1 87 79 T able 3: 204 ST-classes for order 27 the 184 generators that are not self-dual, we obtain 388 generators for the 12911 H-classes. Details of the 204 kno wn ST-classes are summarised in T able 3. Twen ty of these ST-classes are self-dual; the remaining 184 each split into t wo S- classes. In the table, the columns headed “ ||C || S T ” giv e the size of an ST- class C , the next columns “#” giv e the num b er of suc h classes, and the columns “#split” giv e the n umber of these that split into t wo S-classes. There is one “gian t” class G of size 5765 HT-classes (11483 H-classes) and 203 small classes (maximum size 36 HT-classes). T am ura’s matrix R generates the third-largest class, of size 33 HT-classes (66 H-classes). Unlike order 26 (see § 4), there is no ob vious sub division of the classes in to t yp es. Automorphism group orders for the 12911 H-classes are summarised in T able 4. T amura’s matrix R has group order 3. Order m ultiplicity 1 12738 2 26 3 131 6 16 T able 4: Distribution of group orders for 12911 H-classes 6 Results for order 33 F or order 33, the Barba b ound (1) gives D (33) < 516 × 2 74 . Using the algorithm describ ed in [3], we hav e sho wn that none of the 13670 candidate Gram matrices G satisfying det( G ) 1 / 2 ≥ 470 × 2 74 can decomp ose into a pro duct RR T , where R ∈ { +1 , − 1 } 33 × 33 . Th us, w e ha ve D (33) < 470 × 2 74 . 11 Figure 2: ST-classes of size 28 and 36 for matrices of order 27 On the other hand, Solomon [21, 18] found a matrix R ∈ { +1 , − 1 } 33 × 33 with det( R ) = 441 × 2 74 , which is greater than 0 . 9382 of the upp er b ound. Th us, w e know that 441 ≤ D (33) / 2 74 < 470 . It is plausible to conjecture that the low er b ound is b est p ossible and D (33) = 441 × 2 74 . Proving this seems difficult, for reasons giv en in [3, § 7.1]. In this section we find a large num b er of H-classes of solutions with determinan t 441 × 2 74 . Ev en if the conjecture pro ves to b e incorrect, the same tec hniques should b e applicable to find many or all H-classes of solutions with larger determinan t. Starting from the Gram matrix G = R T R = RR T corresp onding to Solomon’s { +1 , − 1 } matrix R , our randomised tree searc h algorithm can find many solutions with the same determinant. Using row/column switching and duality , we can find a huge n umber of inequiv alent solutions. F or example, starting from Solonon’s matrix R and iterating the op eration of row switc hing only , we found 37030740 H-classes in 11 iterations b efore stopping our program b ecause it was using to o muc h memory . Clearly a different strategy is needed. 12 Exploring the switc hing graph using random walks Giv en solutions A 0 and B 0 , we can generate random w alks ( A 0 , A 1 , A 2 , . . . ) and ( B 0 , B 1 , B 2 , . . . ) in the graph defined b y ro w/column switching and transp osition. Each vertex on a walk is connected by a row/column switch- ing op eration and p ossibly transp osition to its successor. If A 0 and B 0 are in different connected comp onen ts, then the tw o ran- dom w alks can not in tersect. Ho wev er, if A 0 and B 0 are in the same con- nected comp onen t, of size s say , then we exp ect the tw o walks to in tersect ev entually , and probably after O ( √ s ) steps unless the mixing time of the w alks is to o long (this dep ends on the geometry of the comp onent, which is unkno wn). Our implemen tation uses self-av oiding random w alks. Eac h walk is stored in a hash table so we can quickly chec k if a new vertex has already b een encountered in the same walk (in which case we try one of its neigh- b ours) or in the other w alk (in whic h case we hav e found an intersection). If, during a walk, all neighbours of the current vertex hav e b een visited, then it is necessary to backtrac k. This o ccurs rarely since the mean degree of a v ertex is large (see be lo w). W e fix A 0 = R and choose B 0 randomly . Usually (ab out 90% of the time) R and B 0 are in the same connected component, nearly alw ays the “gian t” comp onent G . Otherwise, B 0 is in a “small” component (of size sa y s ) and we disco ver this b y being unable to contin ue the self-av oiding w alk from B 0 past B s − 1 . In this w a y we find man y vertices of the gian t comp onen t G , and also man y “small” ST-classes. W e can gather statistics from the random walks. F or example, w e would lik e to estimate ||G || S T , the total num b er of HT-classes (i.e. connected com- p onen ts) in the graph, the n umber of ST-classes, the mean degree of eac h v ertex, etc. If implemen ted as describ ed abov e, the random w alks are not uniform o ver the v ertices of the connected comp onen ts containing their starting p oin ts. They are approximately uniform ov er e dges , so the probability of hitting a v ertex v dep ends on the degree deg( v ). W e can either take this in to accoun t when gathering statistics, or av oid the problem b y accepting a candidate v ertex v with probabilit y 1 / deg( v ). In this wa y the vertices are sampled uniformly if the walks are long enough. The dra wback is that we ha ve to compute the degrees of all candidate v ertices (all neighbours of the current v ertex), whic h migh t b e time-consuming. 13 Results from random w alks W e estimate that the ov erall size of the graph is (3 . 08 ± 0 . 09) × 10 9 when measured in HT-classes. (In terms of H-classes the num b ers are roughly doubled, since self-dual classes are rare.) The giant comp onen t G has size ||G || S T ≈ (2 . 83 ± 0 . 08) × 10 9 . In G the mean degree of eac h vertex is ab out 20, so there are ab out 2 . 83 × 10 10 edges. W e estimate that there are ab out 5 × 10 7 small ST-classes, with mean size ab out 5. Of these we found more than 8 × 10 4 so far, with the largest ha ving size 2136 (see T able 5). W e hav e found 1639 singletons (ST-classes of size 1) and 5412 self-dual matrices. One self-dual singleton w as found. The automorphism group orders observed during random w alks are in { 1 , 2 , 4 } , with orders greater than 1 b eing rare. Solomon’s R has group order 2. Although most of these observ ations are imprecise, since they dep end on random sampling, w e can at least claim the follo wing: Theorem 3. F or or der 33 and determinant 441 × 2 74 , the ST-class G gen- er ate d by Solomon ’s matrix R is self-dual and has size ||G || S T > 197 × 10 6 . Ther e ar e at le ast 8 × 10 4 smal ler ST-classes, 20 of which ar e liste d in T a- ble 5 . Pr o of. As usual, the proof is computational. Starting from R , we found 197122852 HT-classes in 9 iterations of row/column switc hing and taking duals. Starting random w alks from R and R T , and p erforming ro w/column switc hing only , we found an in tersection. Thus R ↔ R T . It follows that G is self-dual. Remarks 1. Although R is not self-dual, we found man y self-dual matrices in the gian t class G in the course of v arious random walks. Eight suc h matrices are at distance 3 from R . The existence of a self-dual matrix in G is sufficient to show that G is self-dual. 2. In addition to the giant class G of size ab out 2 . 83 × 10 9 , we found 20 ST-classes of size ≥ 900, as listed in T able 5, with the largest ha ving size 2136. Classes of the same size marked “(a)” and “(b)” are different, so there are (at least) tw o disjoin t classes of size 999. 3. T able 5 giv es the num b er of times λ i that we found the same class of size s i . This statistic is mentioned b ecause it indicates how well we ha ve 14 size s i λ i size s i λ i 2136 2 1100 1 1300 4 1069 2 1276 2 1011 1 1246 1 1008 1 1205 4 999 1(a) 1188 4 999 2(b) 1187 1 993 2 1148 2 958 2 1134 2 918 3 1104 2 909 3 T able 5: Some large ST-classes for order 33 sampled the search space. Excluding the giant class G , the 20 largest known classes, of total size P s i = 22888, were found P λ i = 42 times. Consider the union U of these classes as a sample from the space, and assume that the space is sampled uniformly . Excluding the 20 “hits” used to select U , there are P ( λ i − 1) = 22 additional “hits” on U . Thus, the fraction ρ of the space sampled is ρ ≈ P ( λ i − 1) / P s i = 22 / 22888 ≈ 1 / 1040. Under our assumption, the probability P of missing a given class of size ≥ 2136 is b ounded b y P ≤ (1 − ρ ) 2136 < 1 / 7 . (5) On the other hand, random sampling hit the gian t class 1 . 04 × 10 6 times, and the estimated size of this class is 2 . 83 × 10 9 , implying that ρ ≈ 1 / 2720, so the estimate (5) on P should b e view ed with caution. The discrepancy b et w een the tw o estimates of ρ ma y b e caused b y nonuniformit y of sampling and/or by an inaccurate estimate of the size of the giant class. 4. It would b e in teresting to know more ab out the graphs asso ciated with “small” ST classes. W e hav e observ ed one graph of size 3 (“ ∨ ”), t w o of size 4 (“ t ” and “  ”), and only one of size 5 (“kite”). Figure 3 shows an example of eac h size in the range 10 , . . . , 19. 5. The reader may ha ve noticed that the graphs displa yed in Figures 1–3 are bipartite (2-colourable), although ne ato did not dra w them in a w ay that mak es this obvious. Computational exp erimen ts hav e sho wn that most, but not all, of the ST classes considered ab o ve hav e bipartite graphs. In particular, the graph of the giant component for order 33 is not bipartite. 15 Ac kno wledgements The author thanks Judy-anne Osb orn for in tro ducing him to the topic and for many stimulating discussions; Will Orric k and Paul Zimmermann for their ongoing collab oration; Brendan McKa y for his graph isomorphism pro- gram nauty [15], which w as used to chec k Hadamard equiv alence; A T&T for the Gr aphviz graph visualization soft ware, in particular ne ato , which w as used to draw the figures; and the Mathematical Sciences Institute (ANU) for computer time on the cluster “orac”. Figure 3: One ST-class of each size 10 , . . . , 19 for order 33 16 References [1] G. Barba, In torno al teorema di Hadamard sui determinanti a v alore massimo, Giorn. Mat. Battaglini 71 (1933), 70–86. [2] R. P . Bren t, The Hadamard maximal determinant problem , http:// maths.anu.edu.au/ ~ brent/maxdet/ [3] R. P . Brent, W. P . Orrick, J. H. Osb orn and P . Zimmermann, Maximal determinan ts and saturated D-optimal designs of orders 19 and 37, submitted. Av ailable from . [4] R. P . Brent and J. H. Osb orn, General low er b ounds on maximal determinan ts of binary matrices, submitted. Av ailable from http: //arxiv.org/abs/1208.1805 . [5] A. E. Brou wer, An infinite series of symmetric designs , Math. Centrum, Amsterdam, Rep ort ZW 202/83 (1983). [6] R. H. F. Denniston, En umeration of symmetric designs (25 , 9 , 3). In A lgebr aic and ge ometric c ombinatorics , V olume 65 of North-Holland Math. Stud., North-Holland, Amsterdam, 1982, 111–127. [7] H. Ehlich, Determinantenabsc h¨ atzungen f ¨ ur bin¨ are Matrizen, Math. Z. 83 (1964), 123–132. [8] H. Ehlic h, Determinantenabsc h¨ atzungen f ¨ ur bin¨ are Matrizen mit N ≡ 3 mo d 4, Math. Z. 84 (1964), 438–447. [9] J. Hadamard, R ´ esolution d’une question relativ e aux d ´ eterminan ts, Bul l. des Sci. Math. 17 (1893), 240–246. [10] N. Ito, J. S. Leon and J. Q. Longy ear, Classification of 3- (24,12,5) designs and 24-dimensional Hadamard matrices, J. Combin. The ory Ser. A 31 (1981), 66–93. [11] H. Kimura, New Hadamard matrix of order 24, Gr aphs Combin. 5 (1989), 235–242. [12] C. Koukouvinos, S. Kounias and J. Seb erry , Supplementary difference sets and optimal designs, Discr ete Math. 88 (1991), 49–58. [13] S. Kounias, C. Koukouvinos, N. Nikolaou and A. Kakos, The non- equiv alent circulan t D-optimal designs for n ≡ 2 mo d 4, n ≤ 54, n = 66, J. Combin. The ory Ser. A 65 (1994), 26–38. 17 [14] B. D. McKa y , Hadamard equiv alence via graph isomorphism, Discr ete Mathematics 27 (1979), 213–214. [15] B. D. McKay , nauty User’s Guide (V ersion 2.4), Dept. of Computer Science, Australian National Universit y , No vem b er 2009. [16] William P . Orrick, Switching operations for Hadamard matrices, SIAM J. Discr ete Math. 22 (2008), 31–50. [17] William P . Orrick, On the enumeration of some D-optimal designs, J. Statist. Plann. Infer enc e 138 (2008) 286–293. abs/math/0511141v2 [18] William P . Orrick, The Hadamard maximal determinant problem , http://www.indiana.edu/ ~ maxdet/ [19] Judy-anne H. Osb orn, The Hadamard Maximal Determinant Prob- lem , Honours Thesis, Universit y of Melb ourne, 2002, 142 pp. http: //wwwmaths.anu.edu.au/ ~ osborn/publications/pubsall.html [20] R. E. A. C. Paley , On orthogonal matrices, J. Math. Phys 12 (1933), 311–320. [21] B. Solomon, P ersonal comm unication to W. Orrick, 18 June 2002. [22] E. Sp ence, Sk ew-Hadamard matrices of the Go ethals-Siedel type, Cana- dian J. Math. 27 (1975), 555–560. [23] H. T amura, Personal comm unication to W. Orric k, 26 August 2005. [24] Ian M. W anless, Cycle switc hes in latin squares, Gr aphs Combin. 20 (2004), 545–570. [25] A. L. Whiteman, A family of D-optimal designs, Ars Combinatoria 30 (1990), 23–26. [26] W. W o jtas, On Hadamard’s inequalit y for the determinan ts of order non-divisible by 4, Col lo q. Math. 12 (1964), 73–83. [27] C. H. Y ang, On designs of maximal (+1 , − 1)-matrices of order n ≡ 2 (mo d 4), Math. Comp. 22 (1968), 174–180. 18