Abstract
We prove the equivalence of two-symbol supersaturated designs (SSDs) with (even) rows, columns, and , where and and resolvable incomplete block designs (RIBDs) whose any two blocks intersect in at most points. Using this equivalence, we formulate the search for two-symbol -optimal and minimax-optimal SSDs with as a search for RIBDs whose blocks intersect accordingly. This allows developing a bit-parallel tabu search (TS) algorithm. The TS algorithm found -optimal and minimax-optimal SSDs achieving the sharpest known lower bound with of sizes , (16, 26), (16, 27), (18, 23), (18, 24), (18, 25), (18, 26), (18, 27), (18, 28), (18, 29), (20, 21), (22, 22), (22, 23), (24, 24), and (24, 25). In each of these cases, no such SSD could previously be found.
1. Introduction
Two-symbol supersaturated designs (SSDs) are two symbol arrays in which the number of rows is less than or equal to the number of columns. Throughout this paper, an SSD refers to a two-symbol SSD. An SSD with rows and columns is represented by an matrix and will be denoted by or simply by . Each entry of is , and the frequencies of and are the same in each column. Moreover, has no two columns such that or .
The -optimality criterion was defined in [1], for comparing two-symbol SSDs. The criterion compares two arrays of the same size by picking the one that minimize where is the th entry of the matrix for a -array . The term in (1) measures the degree of nonorthogonality between the th and th columns. An SSD is called -optimal if no SSD with the same number of rows and columns having a smaller value exists. For arrays, let and be the frequency of in . The minimax criterion proposed in [1] minimizes first and second. An SSD is called minimax-optimal if no other SSD with the same size has a lower or the same with a smaller (see [2]). In [2], a search algorithm generalizing the exchange algorithm of [3] was provided to construct -optimal and minimax-optimal SSDs. An -optimal and minimax-optimal SSD with 16 rows and 60 columns was found by [4]. In [5], a simulated annealing algorithm was used for finding -optimal and minimax-optimal cyclic SSDs. In [6], tabu search (TS) was used to construct -optimal SSDs with good properties by constructing supplementary difference sets. In [7], a TS procedure for constructing -optimal and minimax-optimal -circulant SSDs was implemented. Recently, all isomorphism classes of -optimal and minimax-optimal -circulant SSDs with rows, columns, and were classified in a computer search by [8]. They also classified all isomorphism classes of -optimal and minimax-optimal -circulant SSDs with (mod 4) and . For a comprehensive review of SSDs, see [9].
In Section 2, we provide some background material on lower bounds and the and minimax optimality of SSDs. In Section 3, we prove an equivalence between SSDs with (even) rows, columns, and , where and and resolvable incomplete block designs (RIBDs) such that any two distinct blocks intersect in at most points. In Section 4, using this equivalence, we formulate the problem of constructing -optimal and minimax-optimal SSDs with as a problem to find RIBDs whose blocks intersect in at most points for . We formulate the construction of such RIBDs as an optimization problem. Unlike many optimization problems where a good approximate solution is sufficient, in the construction of such resolvable designs (as in the construction of other combinatorial designs), the main goal is to find an optimum solution. For this purpose, an algorithm based on TS [10] is developed in Section 5. Sets (blocks) with elements from a set of cardinality are stored as bit strings where the number of bits is equal to . Each bit corresponds to exactly one element of . Thus, a set is represented by a bit string in which the bits corresponding to the elements of that set are and all others bits are . (For example, in , bit string represents the set .) Our data structure for storing sets (blocks) is the same as that in [11]. This allows us to exploit bit-parallelism as in [11] for computing the intersections of the sets (blocks) by using bitwise operations. Thus, all our computations are made using bit-parallel Boolean instructions, which in praxis (on a x86-64 CPU) implies that 64 bits of data are processed at once. This improves the overall performance by a factor of (as is less than 64 in all the cases we studied). In Section 5, we also provide the computational complexity analysis of our algorithm. The implementation details of our algorithm are discussed in Section 6. We end the paper with concluding remarks in Section 7.
The bit-parallel TS algorithm was able to construct fifteen previously unknown -optimal and minimax-optimal SSDs of sizes , (16, 26), (16, 27), (18, 23), (18, 24), (18, 25), (18, 26), (18, 27), (18, 28), (18, 29), (20, 21), (22, 22), (22, 23), (24, 24), and (24, 25). All these SSDs, their values, s, and s, as well as their Gram matrices, are provided in Tables 1–30. The newly found -optimal and minimax-optimal SSDs could not have been found by the algorithms in [2] for despite running these algorithms for a very long time. So, the TS algorithm in this paper outperforms the algorithms at least for the SSD cases searched in this paper. The bit-parallel TS algorithm also found all -optimal and minimax-optimal SSDs obtained by the algorithms in [2]. We made a comparison of TS with the algorithms. This is because in TS, , and , the objective function was chosen with respect to both the and the minimax criteria. In addition, algorithms were successful in locating all -optimal SSDs achieving the Ryan and Bulutoglu lower bound [2] for except the 14-row, 16-column case, where the Ryan and Bulutoglu lower bound for this case was recently improved [12]. No other previous algorithm had solved so many cases.
SSDs are used in computer experiments; software testing; medical, industrial, and engineering experiments; chromatography (separation science); and in biometric applications. In [3], an SSD with runs (rows) and factors (columns) for a crash test experiment on a planned new four-wheel drive range was discussed, where the objective was to find the best possible subset of safety features among the 54 proposed. In [13], an SSD with runs (rows) and factors (columns) to fine-tune 16 potential factors affecting the thermal performance of project homes was proposed, where two additional columns were used as blocking factors. For designing a multistage axial compressor (turbine engine), the design engineer has selected 27 potentially important factors [14]. In [14], SSDs with factors (columns) and , , and runs (rows) from [3] and [15], respectively, were compared by using data from a computer experiment for the design of a multistage axial compressor. Our newly found run (row), factor (column) -optimal, and minimax-optimal SSD could have also been used in this study. In [16], an SSD with runs (rows) and factors (columns) was used in composite sampling for monitoring pesticide residues in water. Our newly found run (row), factor (column) -optimal, and minimax-optimal SSD could also have been used for the same purpose. In [17], an run (row), factor (column) SSD was proposed in place of a Plackett-Burman design that had been previously used in [18] for developing an epoxide adhesive system for bonding a polyester cord. We propose our newly found run (row), factor (column) -optimal, and minimax-optimal SSD for the same purpose.
2. Lower Bounds for and and Minimax Optimality of SSDs
In [3, 19], it was independently shown that
Bound (2) can be achieved only if and (mod 4), or if and (mod 4) for some positive integer . Bound (2) was independently improved in [20, 21]. These improved lower bounds are equal when they both apply.
Let and , where and are the floor and ceiling functions. The following theorem by [2] provides an improved version of the [20] lower bound.
Theorem 1 (see [2]). Let be a positive integer such that . Then, there exists a unique nonnegative integer (which depends on and ) such that and (mod 4). Define . (a)If (mod 4), then(b)If (mod 4), thenwhere for even and for odd
In this paper, we search for -optimal SSDs achieving the lower bound defined by (3), (4), (5), and (6) in Theorem 1 with . The following theorem shows that -optimality is a sufficient condition for minimax optimality if .
Theorem 2 (see [2]). Let be an -optimal SSD. (a)If (mod 4) and , then is minimax-optimal(b)If (mod 4) and , then is minimax-optimal
3. The Equivalence between SSDs and RIBDs
An incomplete block design (IBD) with parameters , denoted by IBD , is a pair where is a -set of points and is a collection of -subsets (blocks) of , . The parameters must satisfy the condition
An IBD with parameters satisfying (7) is called a resolvable incomplete block design, denoted by RIBD, if the collection of blocks can be partitioned into subsets called parallel classes of size , each of which partitions the point set.
Henceforth, will denote a positive even integer greater than or equal to 8. Let and be two different parallel classes on points. `1Define their parallel class intersection matrix (PCIM) as the matrix with entries defined by (see [22]). An IBD has parallel classes. Thus, for an arbitrary fixed parallel class , there are PCIMs of the form . Since each point belongs to exactly one block of a parallel class, both the column and row sums in a PCIM are . Then, by relabeling the blocks in parallel classes if necessary, each of the PCIMs associated with can be assumed to be one of for . For any matrix , let ; hence, .
Theorem 3. An RIBD with parameters such that for any two distinct parallel classes and , exists if and only if a row, column, -array with each column orthogonal to the all s column and , where exists.
Proof. Suppose that are the parallel classes of the RIBD . Then, for each parallel class (), we define the column vector as follows:
for each , where is the th entry of . Note that the frequencies of and are the same in each column constructed from the RIBD. Hence, these columns form a -array with each column orthogonal to the all 1s column. For any two columns and of defined by the parallel classes and , we have by (8) and (9). This implies that , where .
Conversely, suppose that is a array with each column orthogonal to the all 1s column. For each column () of the array , there exist two blocks and that partition , such that if the th entry of is , then is contained in the block ; otherwise, is contained in . Clearly, these two blocks form a parallel class. Since the frequencies of and in each column are both , each block has size . Now, , where implies that for any two parallel classes and defined by the th and th columns.
Next, we provide an example for Theorem 3, where is an RIBD corresponding to a 6 row, 8 column, array with each column orthogonal to the all 1s column and .
Example 1. , , , , , , , , , , , , , , , and .
4. The Optimization Problem
In this section, using the equivalence given in Theorem 3, we formulate the problem of constructing an -optimal and minimax-optimal SSD with (even) rows, columns, and for (i.e., ) as a discrete optimization problem. The following theorem is used to define the objective function for this optimization problem.
Theorem 4. Let be an SSD with rows and columns and be the parallel classes of the RIBD defined by the columns of according to Theorem 3. Let , , and ; then, the following hold: (a)(b)(c), where
Proof. By the proof of Theorem 3 and (8), we have (or ) for some . Then, . Statements (b) and (c) follow from (a).
By Theorem 3, a feasible solution to our optimization problem is a RIBD with parameters . However, since each parallel class of the RIBD is uniquely determined by one of its blocks, a feasible solution reduces to a set of blocks, each of size . Then, based on Theorem 4, we define the objective function as
Computing the intersections of blocks is the bottleneck in computing . We used the bit array data structure to store each set (block) as in [11]. This allowed speeding up the calculation of intersections of blocks by using bit-parallel Boolean instructions. To determine the number of points of intersection of two blocks, we used the SSE4.2 SIMD instruction, _mm_popcnt_u64, included in the recent general-purpose processors. It counts the number of bits set to 1 in a word of 64 bits. Thus, in the C language, is calculated by , where is the binary representation of the block with . These instructions increase the speed by a factor of .
To construct -optimal and minimax-optimal SSDs based on Theorem 2, it is necessary to require that for (i.e., ). For this purpose, we define for , where is the lower bound given in Theorem 1. Then, we modify the objective function (11) to
It follows from Theorem 4 (c) and (12) that if the objective function (13) reaches the value , then we have for . Hence, Theorems 2, 3, and 4 (a) and (b) imply that an -optimal and minimax-optimal SSD with (even) rows, columns, and for is found whenever the objective function (13) achieves for a RIBD.
5. Tabu Search for SSDs
The TS algorithm introduced by [10] is an iterative metaheuristic technique used to search for a solution that minimizes an objective function over a set of feasible solutions . TS has been used successfully to construct -optimal designs, constant weight codes, 1-rotational resolvable balanced incomplete block designs, covering designs, and -optimal and minimax-optimal -circulant SSDs (see [7, 23–26]).
TS is based on a neighborhood search (NS). In NS, each feasible solution has an associated set of neighbors, , called the neighborhood of . It starts with a given initial feasible solution and searches the set by moving from one solution to another in its neighborhood. At each iteration, a move from the current solution to a best one in regardless of whether is made. If more than one solution has the same minimum value, the tie is broken randomly. However, the main shortcoming of NS is cycling through a set of solutions, i.e., keeping on revisiting the same set of solutions. To prevent cycling, TS maintains a list called the tabu list of length. Each move in is removed after iterations.
Sometimes, the tabu list may forbid certain desirable moves, such as those that lead to a better solution than the best one found so far. An aspiration criterion is introduced to cancel the tabu status of a move when this move is judged useful.
TS stops when the objective function reaches the lower bound . However, there is no guarantee of reaching the lower bound, and the search process is stopped if the number of iterations used without improving the best solution exceeds a preset limit.
Two IBDs with parameters are defined as neighbors if they are identical for every parallel class but one, and in that parallel class, there are exactly two points that switch blocks. A swap move is entirely determined by the vector , where points and are switched in the parallel class . The definitions of the neighborhood and the objective function in Section 4 allow calculating the change in the objective function (13) value without recomputing the objective function (13). The only blocks that change after the move are and .
Whenever the points and switch blocks in the parallel class , the tabu list forbids the exchange of the points and at the parallel class in the subsequent iterations. Formally, the tabu list consists of vectors , where the points and were forbidden to be exchanged during the preceding iterations, in the parallel class . The tabu list length was adjusted experimentally. For the problem instances in this paper, the best seems to be some integer between 6 and 8. The pseudocode of our TS algorithm is presented in Algorithm 1.
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||
In the above described algorithm, the computation time is mainly spent on iterations. Hence, we next provide the complexity analysis of each iteration. Let be an RIBD (a feasible solution to our optimization problem). In Algorithm 1, for each block that changes to , there are blocks, that do not change. Let . Then, the objective function (13) is updated according to
Since the intersection of two blocks is performed in bitwise operations, the complexity to update the objective function after a move is . The only blocks that change after the move are and (), and there are possible changes to . Hence, the size of the neighborhood of any RIBD is . Then, the overall time spent for each iteration of this algorithm is
Let denote the expected number of iterations of the algorithm for the column and row case. Then, for , the expected time complexity for each run of this algorithm is . For the most difficult cases, the algorithm was run for 4,000,000 times. So, the overall expected running time of the algorithm was . If we had not used bit-parallelism, then the expected time complexity for each run would have been .
6. Implementation Details
The TS algorithm described above was programmed in C, and all computations were carried out on a 2.67 GHz or 2.4 GHz processor. The source code of the algorithm can be requested by sending an email to the first author.
The TS algorithm was used to construct fifteen -optimal and minimax-optimal SSDs achieving the lower bound of [2] with for and . The existence question for each of these SSDs was previously unknown.
Initial computational experiments showed that the strategy of running the algorithm a larger number of times with a smaller was better than running the algorithm a smaller number of times with a larger . For example, for and , 8,000 runs of the algorithm with found 3 optimum solutions, whereas only 1 optimum solution was found by 2,000 runs with . Then, the TS procedure was carried out at most 80,000 times using on each instance tested. With these values, the TS algorithm did not produce any optimum solutions for . However, the best SSDs found for these cases had for and an value equal to and , where is larger by one than the of a hypothetical -optimal and minimax-optimal SSD achieving the lower bound of [2] with . Then, for these cases, the TS algorithm was carried out at most 4,000,000 times with . Since the TS algorithm runs are completely independent and no information is exchanged, an independent-thread parallelization strategy was used. Different random number seed values were used to avoid an overlapping search.
The bit-parallel TS was able to construct fifteen previously unknown -optimal and minimax-optimal SSDs. Table 31 lists all cases where the best obtained -optimal SSD is minimax-optimal. Each row of this table corresponds to an SSD. The first column shows the number of rows and columns for the SSDs. Column nRUNs gives the number of runs it took the TS algorithm to find an -optimal and minimax-optimal SSD. The last column gives the CPU time. The row, column case was solved in 585.4 CPU hours. However, with the independent-thread parallelization approach, this search took only 16.1 hours, using a cluster with 40 threads.
In Table 31, we do not observe that the CPU time always increases with the number of columns or rows. There are also big increases on the CPU times, just by the addition of one column. There are two reasons for these observations. Firstly, the geometry of the problem can change as the number of columns of the sought after SSD increases. In particular, let and be the ratio of the number of all row, column, -optimal and minimax-optimal SSDs to the number of all and all locally optimum row, column SSDs. It is possible that for fixed , and/or is not a nonincreasing function of . This is mainly because not every -optimal and minimax-optimal, column, row SSD can be obtained by adding a column to an -optimal and minimax-optimal, column, row SSD. Then, the probability of finding an -optimal and minimax-optimal SSD in one iteration of the TS algorithm may actually increase going from columns to columns. Secondly, TS is not a deterministic algorithm and the CPU times are random with potentially large variances. Large variances may easily blur an increasing pattern. In fact, this is more of an issue for the previously unsolved difficult cases.
The only other algorithm that is competitive with the TS algorithm is the algorithm. For the algorithm, each new random starting SSD is independently picked from the previous random starting SSDs. So, each trial of the algorithm with a new random starting SSD can be thought as a Bernoulli trial with a success probability of of finding an SSD that achieves the best known and minimax lower bounds. Then, we can use in [2] as a surrogate for the time it takes to find an SSD that achieves the best known and minimax lower bounds. Now, where
For the most difficult cases, is very small, making very large. So, it is possible to get very lucky and find a solution very quickly or be very unlucky and not be able to find a solution after a very long time. Hence, for the most difficult cases, we do not gain as much information by comparing CPU times as in the case of exact algorithms. When TS is used, we do not know the distribution of the random variable as the independence of trials is no longer a valid assumption. Finding an empirical distribution for would require repeating the computational experiments in the paper a large number of times. This is neither feasible due to resource and time constraints nor worthed as all cases solved by the TS algorithm have no corresponding CPU times based on the algorithm (or any other algorithm). No corresponding CPU times exist because the algorithm failed to solve them despite being run for a very long time.
7. Concluding Remarks
In this paper, we developed a heuristic algorithm for finding -optimal and minimax-optimal SSDs that is more effective than the previously known most effective algorithm for the same purpose. Our algorithm brings fifteen cases of -optimal and minimax-optimal SSDs within computational reach by taking advantage of the equivalence between an SSD and an RIBD described in the proof of Theorem 3 and bit-parallelism from [11].
Data Availability
All the data is in the paper.
Disclosure
The views expressed in this article are those of the authors and do not reflect the official policy or position of the United States Air Force, Department of Defense, or the US Government. This paper is published on Arxiv [27].
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors thank the High Performance Computing Lab at IIMAS-UNAM for providing computing resources.