Using simulated annealing for locating array construction

Context: Combinatorial interaction testing is known to be an efficient testing strategy for computing and information systems. Locating arrays are mathematical objects that are useful for this testing strategy, as they can be used as a test suite that enables fault localization as well as fault detection. In this application, each row of an array is used as an individual test. Objective: This paper proposes an algorithm for constructing locating arrays with a small number of rows. Testing cost increases as the number of tests increases; thus the problem of finding locating arrays of small sizes is of practical importance. Method: The proposed algorithm uses simulation annealing, a meta-heuristic algorithm, to find locating array of a given size. The whole algorithm repeatedly executes the simulated annealing algorithm by dynamically varying the input array size. Results: Experimental results show 1) that the proposed algorithm is able to construct locating arrays for problem instances of large sizes and 2) that, for problem instances for which nontrivial locating arrays are known, the algorithm is often able to generate locating arrays that are smaller than or at least equal to the known arrays. Conclusion: Based on the results, it is concluded that the proposed algorithm can produce small locating arrays and scale to practical problems.


Introduction
This paper proposes to use simulated annealing to generate locating arrays [1]. A locating array can be used in combinatorial interaction testing, a well-known strategy for testing computing and information systems [2]. Combinatorial interaction testing aims to test all value combinations among some subsets of factors. The most common form of combinatorial interaction testing is t-way testing, which tests all value combinations among every set of t factors. The value of t is referred to as strength and a value combination among t factors is said to be a t-way interaction or an interaction of strength t.
As a running example, let us consider a simple model of a printer (Table 1). The model consists of four factors: layout, size, color, and display. The first three factors take two values, whereas the last factor takes three values. For example, the layout factor takes either portrait or landscape. A test for this model is a four-tuple where each element is a value of each factor; thus there are a total of 24 possible tests.
Covering arrays are mathematical objects that play a key role in t-way testing. A covering array of strength t can be regarded as a set of tests such that every t-way interaction appear at least once. Figure 1 shows a two-way covering array for this printer example. Each row represents a test. From this figure, it can be seen that every combination of values for two factors occurs in at least one of the six rows. For example, among the size and duplex factors, a total of six interactions of strength two are possible, namely, (A4, oneside), (A4, shortedge), (A4, longedge), (A5, oneside), (A5, shortedge), and (A5, longedge). Each one of these interactions appears somewhere in the array. Testing cost increases as the number of tests increases; thus the problem of finding a small covering arrays is of practical importance. Indeed, there has been a large body of research on generation of covering arrays of small size.  Locating arrays, which were proposed by Colbourn and McClary [1], add more values to covering arrays. For example, a (d, t)-locating array enables to locate any set of interaction faults if there are at most d faults and all the d faults are of strength t [1]. Figure 1 shows an example of a (1, 2)-locating array for the printer model. Now suppose that the outcome of executing a test is either pass or fail and that the outcome is fail if and only if an interaction fault appears in the test. The array enables to locate any fault if it is a two-way interaction fault (t = 2) and there is no other fault (d = 1). For example, if the set of failing tests is {4, 5, 10}, then the faulty interaction is (A5, No), because no other two-way interaction yields this set of failing tests.
In spite of the useful ability described above, studies on locating arrays are still in an early stage and there are few ways of constructing locating arrays. The paper proposes a simulated annealing-based algorithm for generating (1, t)-locating arrays. The algorithm uses this meta-heuristic algorithm for finding a locating array of a fixed size. By repeatedly running simulated annealing with different array sizes, the algorithm constructs a locating array of a small size. We conducted a set of experiments to evaluate the performance of the proposed algorithm. The results obtained show that the algorithm scales to problems of practical sizes when strength t is two and that it is able to produce locating arrays of optimal or suboptimal sizes for small problem instances.
This rest of the paper is organized as follows. Section 2 presents the definition of locating arrays. Section 3 shows the simulated annealing algorithm for finding a locating array of a given size. Section 4 describes the whole algorithm which uses the simulated annealing algorithm as a subroutine. Section 5 describes a set of experiments we have conducted and their results. Section 6 summarizes related work. Section 7 concludes the paper.

Locating arrays
The System Under Test (SUT) is modeled as where v i represents the number of values that can be taken by the ith factor. Unless otherwise explicitly mentioned, we mean by an array A = [a ij ] a matrix that has k columns and satisfies the condition a ij ∈ {0, ..., v j − 1}. A matrix of no rows is also an array.
An interaction T is a set of factor-value pairs each having a different factor. Formally, an interaction is a possibly empty subset of {(j, V j ) : j ∈ {1, ..., k}, V j ∈ {0, ..., v j −1}} such that no two distinct elements (j, V ), (j ′ , V ′ ) ∈ T have the same factor, i.e., j = j ′ . An interaction T is t-way or of strength t if and only if |T | = t (0 ≤ t ≤ k). An interaction of strength 0 is an empty set. We say that a row (b 1 , b 2 , . . . , b k ) of an array covers an interaction T if and only if b j = V for any (j, V ) ∈ T .
Given an array A, we let ρ A (T ) denote the set of rows that cover interaction T . Similarly, for a set T of interactions, we let ρ A (T ) = T ∈T ρ A (T ). Also we let I t denote the set of interactions of strength t. For example, suppose that the SUT model is (2, 2, 2), which represents a system with three factors that have two values. Then I 2 contains, for example, {(1, 0), (2, 1)}, Colbourn and McClary introduce several versions of locating arrays. Of them the following two are of our interest: Definition 2 An array A is (d, t)-locating if and only if In these definitions, ρ A (T 1 ) = ρ A (T 2 ) ⇔ T 1 = T 2 is equivalent to trivially holds for any T 1 , T 2 by the definition of ρ A (·). Hence, what effectively matters is only the other direction; i.e., The intuition behind this definition is as follows. Let us assume that an interaction is either failure-triggering or not and that the outcome of the execution of a test is fail if the test covers at least one failure-triggering interaction; pass otherwise. Given the test outcomes of all tests, a (d, t)locating array or (d, t)-locating array enables to identify any d or any at most d fault-triggering interactions of strength t.

Example
The following matrix is a mathematical representation of the (1, 2)-locating array shown in Figure 1. The SUT is modeled as In the practice of software and systems testing, we are usually more interested in the case where the number of fault-triggering interactions is at most one than in the case where the number of faults is exactly one or more than one. For this practical reason, we restrict ourselves to constructing (1, t)-locating arrays.
Note that I t = ∅ if |I t | = 0. Hence the necessary and sufficient condition for an array A to be (1, t)-locating can be divided into two parts: one for the case where T 1 = ∅ and |T 2 | = 1 and the other one for the case where |T 1 | = |T 2 | = 1. That is, an array A is (1, t)-locating if and only if the two conditions hold: 3 Simulated annealing 3

.1 Overview
Simulated annealing is a class of meta-heuristic search algorithms. In our approach, simulated annealing is used to find a locating array of a given size. A small locating array is generated by repeatedly using simulated annealing while varying the array size. How to vary the size will be discussed in Section 4.
The pseudo-code in Figure 3 shows the overview of our simulated annealing algorithm. The input of the algorithm is the array size, i.e., the number of rows m, and is given as the formal parameter of the function. The number of possible values for the factors, i.e., (v 1 , . . . , v k ), and strength t are not explicitly described as the input, since they never change in the course of repetition of the algorithm. The output is either a (1, t)-locating array obtained or ⊥ in which case the algorithm failed to find a locating array of size m.
The algorithm starts with an m × k array with all entries randomly selected. The temperature T , which is a real value used for directing state search, is set to the initial value t init . Then, an evolution step is repeated until a locating array is found or the number of repetitions reaches a predefined limit k max . An evolution step begins with a random selection of another array A ′ from the set of neighbors of the current array A. A neighbor of an array A is an array that is slightly different from A. The definition of neighbors and the way of selecting a neighbor can make significant effects on search performance.
Cost function f (·) evaluates A ′ with respect to how far it is from a locating array. Our proposed cost function, which will be described in Section 3.2, always takes a non-negative real value and reaches 0 when A is (1, t)-locating. The cost function directs the search, as it is used to determine whether or not to perform the change of the current array. If the change from A to A ′ does not increase the cost, i.e., ∆ = f (A ′ ) − f (A) ≤ 0, then the change is always accepted. Otherwise, the change is performed with probability exp −∆/T , where T is the temperature. This probabilistic decision is intended to avoid trapping in local minima. Then, the temperature is decreased and the evolution step is repeated. Initialize A ⊲ A is an m × k array randomly generated 5: T ← t init ⊲ T represents the temparature.

Cost function
The cost function must be such that its value decreases if the current solution (array) becomes closer to a locating array. As stated in Section 2, the necessary and sufficient condition that an array is (1, t)-locating consists of two parts. Based on this, we design a cost function as a sum of two functions, each evaluating each of the two parts. The cost function f () is of the form: where f 1 (A) and f 2 (A) correspond to Conditions 1 and 2, respectively, and weight is a control parameter which takes a non-negative real value.
We define: In words, Randomly select an entry a i,j from A 5: We define f 2 (T ) as the number of t-way interactions that share identical covering rows with other interactions. Formally,

Neighbors and neighbor selection
We consider two different strategies for selecting a neighbor of the current solution, i.e., array A. Each of the strategies depends on each definition of neighbors. One is straightforward (Figure 4). We regard A ′ = [a ′ ij ] as a neighbor of A = [a ij ] if and only if A ′ is identical to A except for one entry; i.e., {(a ′ i,j , a i,j ) : a ′ i,j = a i,j } = 1. One neighbor is selected as the next solution from all such neighbors uniformly randomly.
The other strategy is intended to lead the search to an optimal solution more directly. In this strategy, an array A ′ is regarded as a neighbor of the current array A only if 1) ρ A ′ (T ) = ∅ for some T ∈ I t such that ρ A (T ) = ∅, . The first condition states that A ′ covers some interaction T that is not covered by A. The second condition signifies that for some two interactions T 1 and T 2 , their covering rows are the same for A but are different for A ′ . The purpose of restricting the candidates for the next array to these neighbors is to increase the likelihood that the next array will have a reduced cost.
The algorithm of this strategy proceeds as follows ( Figure 5): First, if there is at least one interaction that is not covered by any row in A, the algorithm randomly picks one of such interactions, say T , and randomly selects one of the rows of A. Then it overwrites the row with the interaction T . if I = ∅ then 6: Randomly select an interaction T from I 7: Randomly select a row i 8: A ′ ← A with row i being overwritten with T 9: else 10: if |ρ A (T )| > 1 and a random Boolean value is true then 12: Randomly select a row i such that i ∈ ρ A (T )

13:
Randomly select a factor j involved in T

14:
A ′ ← A with a i,j being altered 15: else 16: Randomly select a row i such that i ∈ ρ A (T )

17:
A ′ ← A with row i being overwritten with T If all interactions of strength t appear in the current solution A, randomly pick a t-way interaction T such that the set of rows where it is covered is identical to the set of rows that cover another t-way interaction T ′ . That is, ρ A (T ) = ρ A (T ′ ) holds. We want to have the next solution A ′ such that To this end, one of the two ways is taken to perform randomly. One is to select a row that covers T and and randomly alter the value on one of the factors of T . The other one is to select a row that does not cover T and overwrite it with T .

Algorithm for constructing locating arrays
In this section, we present the whole algorithm for constructing locating arrays. This algorithm repeatedly executes simulated annealing presented in the previous section with array size m being varied.
Binary search fits for the purpose of systematically changing the value of m. The input of the binary search is lower and upper bounds on the size of the smallest locating array, respectively. We let low and high denote the bounds. The first run of simulated annealing is executed for a size ⌊(low + high)/2⌋. while low ≤ high do 6: size ← ⌊(low + high)/2⌋

7:
A ′ ← SA(size) 8: if A ′ =⊥ then If a locating array is found, then high is updated to ⌊(low + high)/2⌋ − 1; otherwise low is updated to ⌊(low + high)/2⌋ + 1. This process is repeated until low > high. Figure 6 shows the binary search algorithm. The algorithm varies the array size and outputs the smallest locating array among the ones that can be obtained during the search.
As simulated annealing is probabilistic, the binary search may fail in finding a locating array even if there is indeed a locating array of size that falls between the given lower and upper bounds. In such a case, the output is ⊥.
Similarly, each run of simulated annealing can fail in obtaining a locating array of a given size. Having these in mind, we propose a two-phase algorithm as shown in Figure 7. The first phase is a fall-back of the binary search. In this phase, binary search is repeated until a locating array is obtained. The second phase is used to improve the solution that has already been obtained. In this phase, simulated annealing is repeated at most max times for the size smaller by one than the array obtained in the first phase. If a locating array is obtained during the max runs, the size is further decreased by one and the phase is repeated. The termination of the algorithm is ensured by timeout. If timeout occurs, the algorithm will output the smallest locating array found during its execution. The output is ⊥ if no locating array has been obtained.
One problem that remains to be solved is how to compute the lower and upper bounds (Line 4 in Figure 7). At this moment, we obtain these bounds using the lower bound formula provided by Tang, Colbourn, and Yin [3].
The formula is: This lower bound formula assumes that the number of values that a factor can take is the same for all factors; i.e., all v i are the same v. Although in our case, v i can be different for different i, this formula can be used to obtain rough bounds as follows. We select the smallest v i , denoted as min{v i }, and then apply that value to the formula to obtain a lower bound. To obtain an upper bound, we select the greatest v i , denoted as max{v i }, and apply max{v i }+1 to the formula. Strictly speaking, the upper bound obtained this way is merely an approximation in the sense that there is no mathematical guarantee that the value is actually an upper bound. In practice, however, it has turned out that our approximation indeed serves as upper bounds for all problem instances tested in the experiments which will be described in the next section.

Experimental results
We developed a program that implements the proposed algorithm using the Java language. This section presents the results of experiments we conducted using the program. All experiments were performed on a Windows 10 machine equipped with 16 GB memory and a 2.0 GHz Xeon E6540 CPU. The program and problem instances used in the experiments are available at Github at https://github.com/tatsuhirotsuchiya/SA4LA.

Experiment 1
We conduct three sets of experiments. The purpose of the first set is twofold. First, it is intended to compare the two neighbor selection strategies of simulated annealing. As stated in Section 3.3, we consider two different approaches. The first one selects an array that is different from the current solution in a single entry as the next solution. The second one chooses the next solution by "overwriting" the current array with an interaction so that the search can fast converge to a locating array. Thereafter, we call the first strategy baseline and the second one the proposed strategy. The second purpose is to determine the values of control parameters for the remaining experiments. While performing the comparison of the two neighbor selection  strategies, we tested a number of different settings of control parameters. Hence the results can also be used for the second purpose.
To this end, we selected one problem instance, called spin-s, from a wellknown benchmark for combinatorial interaction testing research [4,5,6]. This problem instance is the smallest among the 35 instances in the benchmark. Note that this choice favors the baseline strategy, because the strategy searches the solution space in finer steps than the proposed one. We considered the case t = 2; that is, we considered the problem of constructing (1, 2)-locating arrays.
We tested a total of 64 parameter settings for each strategy. Specifically we tested weight ∈ {0.5, 1.0, 2.0, 4.0}, t init ∈ {0.5, 1.0, 2.0, 4.0}, k max ∈ {2048, 4096, 8192, 16384}. The values of r (the temperature cooling rate in Figure 3) and max (the maximum number of repetitions of Phase 2 in Figure 7) were set as follows: r = 0.999, and max = 3. For each of the settings, we ran the program five times and measured the average size of the resulting locating arrays and the average computation time. Figure 8 summarizes the  results. Different dot types represent different strategies. From the results, it can be clearly seen that the proposed strategy outperformed the baseline strategy. Therefore, we only consider the proposed strategy in the remaining sets of experiments. The circles represent the results of the proposed strategy. They were widely distributed horizontally but fairly concentrated vertically. This means that computation time, which highly depends on the value of k max , does not make significant contribution to reducing the array size. Hence for the remaining part of experiments, we set k max to 2048, the smallest value of the four values we tested. Other parameter values were taken from the case when the smallest size had been obtained when k max = 2048; namely, weight = 4.0 and t init = 0.5.

Experiment 2
The purpose of this set of experiments is to examine if the algorithm scales to large problems that arise in practice. To this end, we applied the proposed algorithm to 35 problem instances presented in [4]. This collection of problem Table 2: Performance in generating (1, 2)-locating arrays. The algorithm was repeated five times for each problem instance. The unit of time is second. instances has been used in several studies in combinatorial interaction testing (e.g., [4]). Table 2 lists these instances. The first column shows the problem name. In the second column, f l 1 1 f l 2 2 . . . represents that the problem instance (SUT model) has l i factors that take f i values (i = 1, 2, . . .).
As in the first set of experiments, we ran the algorithm for each instance five times. We set timeout period to one hour for each run. Compared to the problem instance used in the first experiment, many of the instances were large in size, so that the algorithm often failed to complete its execution within the timeout period. Nevertheless, it was often the case that a locating array was obtained, in which case at least one locating array was constructed during algorithm execution.
Having this in mind, we summarize the results obtained in Table 2 for the case t = 2 and Table 3 for the case t = 3. Each row corresponds to each of the 35 problem instances. For each problem instance, we show two numbers, x and y, in the form of x/y/5. Here x is the number of runs that were completed within the timeout period, while y is the number of runs in which at least one (1, t)-locating array was obtained. Hence 0 ≤ x ≤ y ≤ 5 always hold.
The forth and five columns from left show the total running time and size of the resulting locating array averaged over the x finished runs. The sixth column labeled "min" shows the size of the smallest array among these x locating arrays.
Similarly, the time spent from the start until the last locating array was found and the size of that array are averaged over the y runs and shown in the third and second columns from right. (Note that the algorithm needs to continue to run after the last locating array obtained.) The rightmost column shows the size of the smallest array obtained in the y runs.
The results shown in Table 2 shows that the proposed algorithm scales well to these large problems when the strength t is two. To our knowledge, no previous studies report locating arrays of strength two that are as large (in terms of the number of columns) as those obtained in this set of experiments.
When strength t is three (Table 3), our program failed to generate a locating array for many of the problem instances, mainly because the program ran out of memory. This seems to suggest that further improvement in algorithm design or implementation is required. Still, the program was able to produce locating arrays for 11 instances. To the best of our knowledge, there are no previous studies that reported locating arrays of strength three whose sizes are comparable to those we obtained in this set of experiments.

Experiment 3
The third experiment campaign aims to show how close the locating arrays obtained by the proposed approach are to the smallest-known arrays with respect to the number of rows. At present, to our knowledge, the concrete sizes of locating arrays are known for only a small number of problem instances. Table 4 lists several of such (1, 2)-locating arrays, where the problem instance, known array size, and source reference are shown for each array in the three left columns. We applied the proposed approach to these problems and obtained the results summarized in the remaining columns in the same table, as in Tables 2 and 3. Since these problems are relatively small, the algorithm was always successful in completing its execution within one hour, except for the largest instance. Even for that instance, the algorithm was able to produce a locating array in every run of the program.
Numbers in bold font represent the smaller array size between the known locating array or the one obtained in the experiments. The known array sizes underlined are proved to be optimal (minimum). Such optimal arrays have been known only for small problems. The results show that the arrays obtained by applying our algorithm to these small problems were also optimal or very close to optimal.
For larger instances, our algorithm was able to find locating arrays that are considerably smaller than known ones. The size reduction achieved is substantial, particularly for the first two instances which have arisen from real-world applications [7,8]. The (1, 2)-locating array of size 421 for the first instance [7] was "crafted from tools for covering arrays, using simple heuristics" [11]. The array we obtained using the single algorithm is smaller than this previously known array by more than 30 percent. For the second instance, the locating array we found is also approximately 30 percent smaller than the previously known one.

Related work
Colbourn and McClary proposed locating arrays in [1]. They suggested a few possible applications of locating arrays, including combinatorial interaction testing, group testing, etc. Locating arrays were used for screening interacting factors in wireless network simulator and testbed in [7,12]. Applications to combintatorial interacting testing were discussed in [13].
Mathematical aspects of locating arrays have been explored in [3,14,15]. These studies assume that the number of values that a parameter can take is the same for all parameters. In [3], Tang, Colbourn, and Yin proved a lower bound on the size of locating arrays. This bound was used in our proposed algorithm, as described in Section 4. In [14], an infinite but sporadic series of minimum (2, t)-locating arrays are provided. In [15], the constructions of the minimum (1, 1)-, (1, 1)-, (1, 1)-, and (1, 1)-locating arrays are presented. Mathematical constructions of (1, 2)-locating arrays are proposed in [11]. Computational constructions can naturally handle the case where the domain sizes of parameters are different for different parameters. In [16], a "one-test-at-a-time" greedy heuristic algorithm is proposed to construct (1, 2)-locating arrays. In our previous work [9,10], we proposed an approach that uses a Constraint Satisfaction Problem (CSP) solver to find locating arrays. The arrays obtained by using the CSP-based approach were compared with those obtained by the simulation annealing-based algorithm in Section 5. Recently, techniques based on resampling [17] have been investigated for locating array construction. Studies in this line include [18,19].
In contrast to the small amount of research on computational construction of locating arrays, many approaches have been applied to the problem of generating covering arrays. These approaches include simulated annealing [4,20], as well as greedy heuristics [21,22], tabu search [5,23], genetic algorithms [24], particle swarm [25,26], and others [27,28]. Our proposed algorithm for generating locating arrays shares some similarities with the simulated annealing-based covering array constructions but is different in many aspects, not only in technical details but also in fundamental algorithm designs. For example, the cost function of simulated annealing and neighbor selection strategies are completely different.

Conclusion
This paper proposed a simulation annealing-based algorithm to construct (1, t)-locating arrays. In the algorithm, simulation annealing is used to find a locating array of a given size. By repeatedly running simulated annealing with the size being systematically varied, the algorithm eventually yields a locating array of a small size. Experimental results showed that when t = 2, the proposed algorithm scales to large-sized problems and is able to produce locating arrays that are substantially smaller than those previously known. The results also showed that for small problems, the (1, 2)-locating arrays produced by the algorithm were minimum or close to minimum.
Future research includes several directions. For example, improving the scalability of the algorithm deserves further study, since as the experimental results shows, the current algorithm is often unable to produce a locating array when the strength t of the array exceeds two. More fine-grained tuning of parameter values of simulated annealing may also require further research. Machine learning techniques might be useful for this problem. Another possible direction is to adapt the proposed algorithm to construction of other related mathematical objects that can be used for fault localization. Examples of such mathematical objects include detecting arrays [1], error locating arrays [29], and constrained locating arrays [30].