-In this paper, we outline MATHCHECK2, a combination of a SAT solver and a computer algebra system (CAS) aimed at finitely verifying or counterexampling mathematical conjectures. Using MATHCHECK2 we verified the Hadamard conjecture from design theory for matrices up to order 144 and many additional orders up to 168. Also, we provide independent verification of the claim that Williamson matrices of order 35 do not exist, and demonstrate for the first time that 35 is the smallest number with this property. In the course of our work, we discovered over 500 Hadamard matrices which were not equivalent to any matrices in the comprehensive MAGMA Hadamard database. The crucial insight behind MATHCHECK2 is that a combination of an efficient search procedure (like those in SAT solvers) with a domain-specific knowledge base (à la CAS) can be a very effective way to verify, counterexample, and learn deeper properties of mathematical conjectures (especially in combinatorics) and the structures they refer to. MATHCHECK2 can be seen as a systematic parallel generator of structures referred to by the conjecture-under-verification C, and these conjectures are typically of the form “for all natural numbers n, some combinatorial object exists”. MATHCHECK2 uses a divideand-conquer approach to parallelize the search, and a CAS to prune away classes whose members are guaranteed to not satisfy the conjecture C. The SAT solver is used to verify whether any of the remaining structures for each number n satisfy C, and in addition learn UNSAT cores in a conflict-driven clauselearning style feedback loop to further prune the search space. Finally, our latest version of MATHCHECK2 (updated from our previously presented version [5]) uses programmatic methods to learn clauses which are generated on-the-fly by domain-specific constraints.
“Brute-brute force has no hope. But clever, inspired brute force is the future.” – Doron Zeilberger1
Many conjectures in combinatorial mathematics are simple
to state but very hard to verify. For example, a conjecture
like the Hadamard [
In recent years, conflict-driven clause-learning (CDCL)
Boolean SAT solvers [
Fortunately the strengths of modern SAT solvers and CAS are complementary, i.e., the domain-specific knowledge of a CAS can be crucially important in cutting down the search space associated with combinatorial conjectures, while at the same time the clever heuristics of SAT solvers, in conjunction with CAS, can efficiently search a wide variety of such spaces.
The success of the SAT and CAS combination has been
demonstrated in papers on MATHCHECK [
In this paper, we present a refinement of our methods and
an improved implementation of MATHCHECK2. Previously,
we observed that many of the generated SAT instances used
a common set of clauses to help prune away the search
space. With more and more of these pruning results included
the number of these clauses increased and could harm the
performance of the SAT solver. To assist the SAT solver,
we chose to encode some constraints using a programmatic
approach as presented by Ganesh et al. [
This paper makes the following contributions:
1) An improvement of our prototype MATHCHECK2, a
combination of CAS and SAT solvers for finitely
verifying or counterexampling math conjectures. We discuss
three techniques in Section IV that dramatically improve
the search capabilities of the basic MATHCHECK2
methodology. We further show how we assist the SAT
solver by performing certain checks programmatically
using a CAS and thereby prune away parts of the search
space significantly. All three techniques can be adapted
for other conjectures, beyond the ones considered in this
paper. https://sites.google.com/site/uwmathcheck/.
2) Description of a novel algorithm for finding Williamson
matrices of a given order (or showing that none exist).
This algorithm makes use of recent theorems about the
properties of compressed sequences [
The architecture of our proposed MATHCHECK2 system is outlined in Figure 1. At its heart is a generator of combinatorial structures, which uses data provided to it by CAS functions to prune the search space and interfaces with SAT solvers to verify the conjecture-in-question. The generator contains functions useful for translating combinatorial conditions into Problem
Generator SAT
UNSAT
SAT instance Domain-specific constraints assums : vars( ) ! fT; Fg
SAT solver result (Solution or UNSAT core)
SAT+CAS clauses which can be read by a SAT solver. It is possible to substitute the SAT solver with an SMT solver to simplify the encoding process, but this resulted in a too large overhead in our computations. The generator is currently optimized to deal with conjectures which concern Hadamard matrices from coding and combinatorial design theory.
Once the class of combinatorial objects has been determined,
the script accepts a parameter n which determines the size
of the object to search for. For example, when searching for
Hadamard matrices, the parameter n denotes the order (i.e., the
number of rows) of the matrix. The generator then queries the
CAS (we chose MAPLE in our calculations) it is interfaced with
for properties that any order n instance of the combinatorial
object in question must satisfy. However, since our generator
is written in Python, many CAS functionality is provided
by certain modules such as NUMPY [
Once the generator determines the space to be searched it splits the space into distinct subspaces in a divide-and-conquer fashion. Once the partitioning of the search space has been completed, the script generates two types of files: 1) A single “master” file in DIMACS2 format which contains the conditions specifying the combinatorial object being searched for. These are encoded as propositional formulae in conjunctive normal form. An assignment to the variables which makes all of them true would give a valid instance of the object being searched for (and a proof that no such assignment exists proves that no instance of the object in question exists). 2) A set of files which contain partial assignments of the variables in the master file. Each file corresponds to exactly one subspace of the search space produced by the generator.
The main advantage of splitting up the problem in such a way is that it easily facilitates parallelization. Using domain 2For more information on this format, please refer to http://www. satcompetition.org/2009/format-benchmarks2009.html specific knowledge, we partition the search space into different classes of the same mathematical structure, and since these classes are independent of each other, a cluster of SAT solvers can search the space for each partition in parallel.
As a new feature, the generator does not include all clauses for sanity checks into the files. We rather perform the most common checks programmatically in a feedback loop inside the SAT solver. In this way, we reduce the number of needed clauses in each file, and we observed a speedup of MATHCHECK2 when searching for Hadamard matrices.
Furthermore, an UNSAT core (generated by SAT solvers
such as MAPLESAT [
Definition 1. A matrix H 2 f 1gn n, n 2 N, is called a Hadamard matrix, if for all i 6= j 2 f1; : : : ; ng, the dot product between row i and row j in H is equal to zero. We call n the order of the Hadamard matrix.
Two Hadamard matrices H1 and H2 are said to be equivalent if H2 can be generated from H1 by applying a sequence of negations/permutations to the rows/columns of H1, i.e., if there exist signed permutation matrices U and V such that U H1 V = H2.
One prominent class of special Hadamard matrices are those generated by so-called Williamson matrices.
Theorem 1 (cf. [
2) A, B, C, and D commute pairwise (i.e., AB = BA,
AC = CA, etc.); 3) A2 + B2 + C2 + D2 = 4nIn, where In is the identity matrix of order n.
details see [
For practical purposes, one considers A, B, C, and D in the Williamson construction to be circulant matrices, i.e., those matrices in which every row is the previous row shifted by one entry to the right (with wrap-around, so that the first entry of each row is the last entry of the previous row). Such matrices are completely defined by their first row [x0; : : : ; xn 1] and always satisfy the commutativity property. If the matrix is also symmetric then we must further have x1 = xn 1, x2 = xn 2, and in general xi = xn i for i = 1, : : : , n 1. Therefore, if a matrix is both symmetric and circulant its first row must be of one of the two forms [x0; x1; x2; : : : ; x(n 1)=2; x(n 1)=2; : : : ; x2; x1] [x0; x1; x2; : : : ; xn=2 1; xn=2; xn=2 1; : : : ; x2; x1] (1) according to if n is odd or even.
form (1), i.e., one which satisfies xi = xn i for i = 1, : : : , n 1.
Williamson matrices are circulant matrices A, B, C, and D which satisfy the conditions of Theorem 1. Since they must be circulant, they are completely defined by their first row. (In light of this, we may simply refer to them as if they were sequences.)
There are three types of operations which, when applied
to the Williamson matrices, produce different but essentially
equivalent matrices. For our purposes, generating just one of the
equivalent matrices will be sufficient, so we impose additional
constraints on the search space to cut down on extraneous
solutions and hence speed up the search.
1. Ordering: Note that the conditions on the Williamson
matrices are symmetric with respect to A, B, C, and D. In
other words, those four matrices can be permuted amongst
themselves and they will still generate a valid Hadamard matrix.
Given this, we enforce the constraint that
jrsum(A)j
jrsum(B)j
jrsum(C)j
jrsum(D)j;
where rsum(X) denotes the sum of the entries of the first (or
any) row of X. Any A, B, C, and D can be permuted so that
this condition holds.
2. Negation: The entries in the sequences defining any of
A, B, C, or D can be negated and the sequences will still
generate a Hadamard matrix. Given this, we do not need to
try both possibilities for the sign of the rowsum of A, B,
C, and D. For example, we can choose to enforce that the
rowsum of each of the generating matrices is nonnegative.
Alternatively, when n is odd we can choose the signs so they
satisfy rsum(X) n (mod 4) for X 2 fA; B; C; Dg. In this
case, a result of Williamson [
Because the search space for Hadamard matrices is so large, it is advantageous to focus on a specific construction method and use known properties of that construction type to filter the search space. We now list the definitions we need to state the main such filtering theorem we used.
Definition 3. The power spectral density of the sequence A is the nonnegative real sequence
PSDA(s) := jDFTA(s)j2 for s 2 Z where DFTA is the discrete Fourier transform of A. Definition 4. The periodic autocorrelation function of the sequence A is the periodic function given by
PAFA(s) :=
In order to check if a matrix H 2 f 1gn n with rows H0, : : : , Hn 1 is Hadamard, it is necessary to verify that Hi Hj = 0 for all 0 i; j < n with i 6= j. In other words, we want to verify that the componentwise product Hi Hj = hi;0 hj;0 hi;1 hj;1 hi;n 1 hj;n 1 a(d) = aj + aj+d +
j
n = dm which satisfy Then we say that the sequence A(d) = [a(0d); a(1d); : : : ; a(dd)1] is the m-compression of A.
Since PSDX (s) is always nonnegative, equation (3) implies that PSDA(d) (s) 4n (and similarly for B, C, D). Therefore if a candidate compressed sequence A(d) satisfies PSDA(d) (s) > 4n for some s 2 Z then we know that the uncompressed sequence A can never be one of the sequences which satisfies the preconditions of Theorem 2.
In the course of our research we discovered the following useful properties of the compressed sequences which arise in our context. These lemmas significantly reduce the number of SAT instances we need to generate.
Lemma 1. If A is a sequence of length n = dm with 1 entries, then the entries ai(d), i 2 f0; : : : ; d 1g, have absolute value at most m and ai(d) m (mod 2).
symmetric.
IV. ENCODING AND SEARCH SPACE PRUNING TECHNIQUES
An attractive property of Hadamard matrices when encoding them in a SAT context is that each of their entries is one of two possible values, namely 1. We choose the encoding that 1 is represented by true and 1 is represented by false. We call this the Boolean value or BV encoding. Under this encoding, the multiplication function of two x, y 2 f 1g becomes the XNOR function in the SAT setting, i.e., BV(x y) = XNOR(BV(x); BV(y)).
The Williamson encoding is very similar to the general encoding but with fewer variables; we merely have the 4 n+1 2 variables a0; a1; : : : ; ad(n 1)=2e; b0; : : : ; bd(n 1)=2e;
c0; : : : ; cd(n 1)=2e; d0; : : : ; dd(n 1)=2e:
we must enforce the conditions Hj ) = 0 for i 6= j rsum(Ai Aj + Bi Bj + Ci Cj + Di Dj ) = 0 for i 6= j: (4) Like in Section IV-A this is done by defining new variables to represent the entries of the componentwise products. Also, note that because of the symmetry and circulant properties most of the conditions to enforce will be identical. In fact, it is only necessary to encode (4) for i = 0 and j = 1, : : : , n2 1 to ensure that such matrices generate a valid Hadamard matrix.
As a special case of compression, consider what happens when d = 1 and m = n. In this case, the compression of A is a sequence with a single entry whose value is Pn 1 k=0 ak = rsum(A). If A, B, C, and D are f 1g-sequences which satisfy the conditions of Theorem 2, then the theorem applied to this m-compression says that PAFA(1) (0) + PAFB(1) (0) + PAFC(1) (0) + PAFD(1) (0) = 4n which simplifies to
rsum(A)2 + rsum(B)2 + rsum(C)2 + rsum(D)2 = 4n; and by Lemma 1 each rowsum must have the same parity as n.
In other words, the rowsums of the sequences A, B, C, and
D decompose 4n into the sum of four perfect squares whose
parity matches the parity of n. Since there are usually only a
few ways of writing 4n as a sum of four perfect squares this
severely limits the number of sequences which could satisfy the
hypotheses of Theorem 2. Furthermore, some computer algebra
systems contain functions for explicitly computing what the
possible decompositions are (e.g., PowersRepresentations
in MATHEMATICA and nsoks by Joe Riel of Maplesoft [
When using this technique it is necessary to encode constraints on the rowsum of the generating matrices, e.g., rsum(A) = 1. This may be simply done by using a binary adder network on the variables a0, : : : , ad(n 1)=2e. We give the variables which appear twice in the first row of A (due to symmetry) a weight of 2 in the binary adder network so that the rowsum is computed correctly.
Because each instance can take a significant amount of time to solve, it is beneficial to divide instances into multiple partitions, each instance encoding a subset of the search space. In our case, we found that an effective splitting method was to split by compressions, i.e., to have each instance contain one possibility of the compressions of A, B, C, and D. To do this, we first need to know all possible compressions of A, B, C, and D. These can be generated by applying Lemmas 1 and 2. For example, when n = 35 and d = 5 there are 28 possible compressions of A with rsum(A) = 1. Of those, only 12 satisfy PSDA(s) 4n for all s 2 Z. The calculation of PSDA(s) was possible to be performed within Python using NUMPY instead of directly querying a CAS. There are also 12 possible compressions for each of B, C, and D with rsum(B) = rsum(C) = 3 and rsum(D) = 11. Thus there are 124 total instances which would need to be generated for this selection of rowsums, however, only 41 of them satisfy the conditions given by Theorem 2.
Furthermore, if n has two nontrivial divisors m and d then we can find all possible m-compressions and d-compressions of A, B, C, and D. In this case, each instance can set both the m-compression and the d-compression of each of A, B, C, and D. Since there are more combinations to check when dealing with two types of compression this causes an increase in the number of instances generated, but each instance has more constraints and a smaller subspace to search through.
After using the divide-and-conquer technique one obtains a collection of instances which are almost identical. For example, the order 40 instances contained 4185 variables and only at most 184 variables differed between instances. Because of this similarity, a short reason why one instance is unsatisfiable may also apply to other instances which are similar. When this is the case, those instances can immediately be pruned away.
MAPLESAT is one SAT solver which supports UNSAT core generation. Provided a master instance and a set of assumptions (variables which are set either true or false), the UNSAT core contains a subset of the assumptions which make the master instance unsatisfiable. Thus, any other instance which sets the variables in the UNSAT core in the same way must also be unsatisfiable.
There are several constraints shared among the different SAT instances which help prune away large regions of the search space. However, the more constraints we include the more the size of each file increases; the SAT solver has more clauses to handle and this causes an increase of its runtime. In the latest version of MATHCHECK2, we oursourced these checks into a DPLL(CAS) style feedback loop which generates clauses on the fly whenever it encounters that the SAT solver is searching in a space where no solution is provably possible.
Specifically, the possible compressions which were not filtered out by the generator were translated into a set of linear constraints which were passed to the SAT+CAS solver. The CAS would then use these constraints to generate learned clauses as the search progressed.
V. VERIFICATION OF THE NONEXISTENCE OF WILLIAMSON
MATRICES OF ORDER 35
We searched for Williamson matrices of order 35 using the techniques described in Section IV with both 5 and 7compression. Despite the exponential growth of possible first rows of the matrices A, B, C, and D, the described pruning results in 21,674 SAT instances of three possible forms, as described in Figure 2. Each instance has subsequently been checked with several SAT solvers, and each one has been discovered to be unsatisfiable. Using MAPLESAT with UNSAT core generation, 19,356 SAT solver calls were necessary to determine that all instances were unsatisfiable.
Our practice was to have people that were not involved in writing the respective code verify its correctness, and to have domain experts verify the application of the theorems used. Furthermore, our confidence of the correctness of our code was strengthened by the successful discovery of Williamsongenerated Hadamard matrices for all the orders 4n with n < 35. These have been determined to be valid Hadamard matrices by the computer algebra system MAGMA. Programmatic 0.01 0.08 0.05 0.28 0.05 1.51 0.75 6.08 1.08 42.21 VI. EXPERIMENTAL RESULTS ON THE HADAMARD
CONJECTURE
We checked all of the Hadamard matrices we computed for equivalence against those in MAGMA’s Hadamard matrix database. In total, our methods generated more than 500 pairwise inequivalent Hadamard matrices which were also not equivalent to any matrices in this database. We submitted these to the MAGMA team and one can download these on our project website (URL in Section I-A).
Experimental Setup and Methodology: The timings were run on the high-performance computing cluster SHARCNET. Specifically, the cluster we used ran CentOS 5.4 and used 64-bit AMD Opteron processors running at 2.2 GHz. Each SAT instance was generated using MATHCHECK2 with the appropriate parameters and the instance was submitted to SHARCNET to solve by running MAPLESAT on a single core (with a timeout of 24 hours).
Figure 3 contains a summary of the performance of our encoding and pruning techniques. The timings are for searching for Williamson matrices of order n with 25 n 35 and for each of the techniques discussed in Section IV. We did not use Techniques 2 and 3 for orders 29 and 31 as they have no nontrivial divisors to perform compression with, but they were otherwise very effective at partitioning the search space in an efficient way. Technique 3 was effective at cutting down the number of instances generated in certain orders. Although the instances pruned tended to be those which would have been quickly solved, this technique would be especially valuable in a situation where few cores are available, as it would allow the overhead of many SAT solver calls to be avoided.
The timings given in Figure 3 refer to the total amount of time used by MAPLESAT across all the jobs run on SHARCNET for each order and technique. The numbers in parentheses in Figure 3 denote how many MAPLESAT calls returned a result and did not time out. The jobs using the base encoding and Technique 1 which did not time out all returned SAT. All of the jobs using Technique 2 completed without timing out and most the instances were found to be UNSAT. The Technique 3 results were the same as the Technique 2 results except with fewer calls to MAPLESAT as some instances could immediately be determined to be UNSAT.
Figure 4 contains the average runtimes for the SAT instances generated by our script using a compression factor of 2 (and therefore only for even orders). For comparision purposes, MAPLESAT was run both with and without the programmatic functionality and all other techniques enabled.
The idea of combining the capabilities of SAT/SMT solvers
and computer algebra systems or domain-specific knowledge
has been examined by various research groups. Recently, Heule
et al. [
We have successfully presented the advantages of utilizing the power of SAT solvers in combination with domain specific knowledge and algorithms provided by computer algebra systems. Our main mathematical problem was the verification of the Hadamard conjecture for some orders by using MATHCHECK2 to search for and discover Williamson matrices. We verified independently, as requested by D. Ðokovic´, that there is no Hadamard Matrix of order 4 35 which is generated by Williamson matrices. Moreover, we discovered more than 500 Hadamard matrices that are not equivalent to any matrix in the MAGMA Hadamard database.
A future direction is to scale to Hadamard matrices of
higher order. For this, we plan to refine the methods (e.g., by
examining other construction types). We also want to analyze
the UNSAT cores generated by Technique 3 to explain their
effectiveness in certain cases, as well as exploring the usage
of incremental SAT solvers [