How to Estimate the Ability of a Metaheuristic Algorithm to Guide Heuristics During Optimization

TL;DR

This paper introduces a methodology to evaluate how effectively a metaheuristic guides heuristics during optimization by comparing it to a naive placebo version, providing insights into the metaheuristic's specific contribution.

Contribution

It proposes a simple, effective approach using distribution comparison and introduces BER measures to quantify the metaheuristic's guiding ability.

Findings

The methodology successfully distinguishes guiding effectiveness in simulated experiments.

BER measures provide practical significance thresholds for performance differences.

Application to Simulated Annealing and SAT problems demonstrates the approach's utility.

Abstract

Metaheuristics are general methods that guide application of concrete heuristic(s) to problems that are too hard to solve using exact algorithms. However, even though a growing body of literature has been devoted to their statistical evaluation, the approaches proposed so far are able to assess only coupled effects of metaheuristics and heuristics. They do not reveal us anything about how efficient the examined metaheuristic is at guiding its subordinate heuristic(s), nor do they provide us information about how much the heuristic component of the combined algorithm contributes to the overall performance. In this paper, we propose a simple yet effective methodology of doing so by deriving a naive, placebo metaheuristic from the one being studied and comparing the distributions of chosen performance metrics for the two methods. We propose three measures of difference between the two distributions. Those measures, which we call BER values (benefit, equivalence, risk) are based on a preselected threshold of practical significance which represents the minimal difference between two performance scores required for them to be considered practically different. We illustrate usefulness of our methodology on the example of Simulated Annealing, Boolean Satisfiability Problem, and the Flip heuristic.

Tables

Table 1. Table 1: Low, medium, and high levels of SA[Flip]’s parameters.

Parameter	Low	Medium	High
$T_0$	$1$	${10}^2$	${10}^3$
$\alpha$	$0.5$	$0.85$	$0.99$
$M$	$1$	$10$	$20$
$MNI$	$10$	$50$	$100$

Table 2. Table 2: Half-distances of the SA[Flip] parameters, used during calibration.

Parameter	$T_0$	$\alpha$	$M$	$MNI$
Half-distance	$10$	$0.04$	$5$	$10$

Table 3. Table 3: Average percentages of unsatisfied clauses ( Y 𝑌 Y ) in Testing Formulae with different numbers of propositional letters n 𝑛 n , for SA[Flip] and ∅ \emptyset [Flip].

$n$	$\varnothing$[Flip]	SA[Flip]
$50$	$0.000000e+00$	$0.000000e+00$
$75$	$5.256542e-05$	$6.154000e-05$
$100$	$2.587233e-04$	$2.383738e-04$
$125$	$5.583783e-04$	$5.847133e-04$

Table 4. Table 4: Average Percentages of Unsatisfied Clauses ( Y 𝑌 Y ) in Testing Formulae for SA[Flip] and ∅ \emptyset [Flip].

Algorithm	Average $Y$
$\varnothing$[Flip]	$0.0002174168$
SA[Flip]	$0.0002211568$

Table 5. Table 5: Success rates of SA[Flip] and ∅ \emptyset [Flip] for Testing Benchmarks with different number of propositional leters, n 𝑛 n .

$n$	$\varnothing$[Flip]	SA[Flip]
$50$	$100.00\%$	$100.00\%$
$75$	$99.33\%$	$99.12\%$
$100$	$96.04\%$	$96.42\%4$
$125$	$89.96\%$	$89.88\%$

Table 6. Table 6: Success rates of SA[Flip] and ∅ \emptyset [Flip].

Algorithm	Average $Y$
$\varnothing$[Flip]	$96.33\%$
SA[Flip]	$96.35\%$

Table 7. Table 7: Empirical BER values for δ = 0.0000 𝛿 0.0000 \delta=0.0000

$n$	$b^{\ast }$	$e^{\ast }$	$r^{\ast }$
$50$	$0.0000$	$1.0000$	$0.0000$
$75$	$0.0051$	$0.9856$	$0.0093$
$100$	$0.0327$	$0.9164$	$0.0509$
$125$	$0.0638$	$0.8349$	$0.1013$
overall	$0.0254$	$0.9342$	$0.0404$

Table 8. Table 8: Empirical BER values for δ = 0.0100 𝛿 0.0100 \delta=0.0100

$n$	$b^{\ast }$	$e^{\ast }$	$r^{\ast }$
$50$	$0.0000$	$1.0000$	$0.0000$
$75$	$0.0004$	$0.9982$	$0.0014$
$100$	$0.0043$	$0.9870$	$0.0087$
$125$	$0.0017$	$0.9935$	$0.0048$
overall	$0.0016$	$0.9947$	$0.0037$

Table 9. Table 9: Empirical BER values for δ = 0.0200 𝛿 0.0200 \delta=0.0200

$n$	$b^{\ast }$	$e^{\ast }$	$r^{\ast }$
$50$	$0.0000$	$1.0000$	$0.0000$
$75$	$0.0000$	$1.0000$	$0.0000$
$100$	$0.0000$	$1.0000$	$0.0000$
$125$	$0.0000$	$1.0000$	$0.0000$
overall	$0.0000$	$1.0000$	$0.0000$

Equations

$$P(Y|\Pi ,S,{\theta }_{\mathcal{M}},{\theta }_{\mathcal{H}})$$

$$P(Y|M,\Pi ,S,{\theta }_M,{\theta }_{\mathcal{H}})$$

$$P(Y|M=\mathcal{M},\Pi ,S,{\theta }_{\mathcal{M}},{\theta }_{\mathcal{H}})$$

$$P(Y|M=\varnothing ,\Pi ,S,{\theta }_{\varnothing },{\theta }_{\mathcal{H}})$$

$$P(Y_{\mathcal{M}[\mathcal{H}]}<Y_{\varnothing [\mathcal{H}]})$$

$$b=P(Y_{\mathcal{M}[\mathcal{H}]}<Y_{\varnothing [\mathcal{H}]}-\delta )$$

$$r=P(Y_{\mathcal{M}[\mathcal{H}]}>Y_{\varnothing [\mathcal{H}]}+\delta )$$

$$e=P(Y_{\varnothing [\mathcal{H}]}-\delta \le Y_{\mathcal{M}[\mathcal{H}]}\le Y_{\varnothing [\mathcal{H}]}+\delta )=1-b-r$$

$$b^{\ast }=\frac{1}{l\times n^2}\sum _{i=1}^l\sum _{j=1}^n\sum _{k=1}^n{I}_{{\mathbf{Y}}_{\mathcal{M}}^{\ast }[i,j]<{\mathbf{Y}}_{\varnothing }^{\ast }[i,k]-\delta }$$

$$r^{\ast }=\frac{1}{l\times n^2}\sum _{i=1}^l\sum _{j=1}^n\sum _{k=1}^n{I}_{{\mathbf{Y}}_{\mathcal{M}}^{\ast }[i,j]>{\mathbf{Y}}_{\varnothing }^{\ast }[i,k]+\delta }$$

$$e^{\ast }=\frac{1}{l\times n^2}\sum _{i=1}^l\sum _{j=1}^n\sum _{k=1}^n{I}_{{\mathbf{Y}}_{\varnothing }^{\ast }[i,k]-\delta <{\mathbf{Y}}_{\mathcal{M}}^{\ast }[i,j]<{\mathbf{Y}}_{\varnothing }^{\ast }[i,k]+\delta }$$

$$\underset{i=1}{\overset{m}{\bigwedge }}\underset{j=1}{\overset{k_i}{\bigvee }}\pm p_{i,j}$$

$$b^{\ast }=\frac{1}{l}\sum _{i=1}^l{I}_{Z_{\mathcal{M}[i]}<Z_{\varnothing }[i]-\delta }$$

Full text

Keywords: Algorithm Analysis, Metaheuristics, Heuristics, Simulated Annealing, Boolean Satisfiability

1 Introduction

Metaheuristics and heuristics are widely accepted optimization tools within operations research community (Caserta & Voß 2010). They are used to approximately, but efficiently, solve problems that are too hard to be solved using exact algorithms (Nesmachnow 2014).

Heuristics are problem-specific techniques that, in general, quickly find good solutions to given problems, although there are no guarantees that those solutions will always be optimal. Heuristics can be used only to solve problems for which they have been specifically designed. Metaheuristics, on the other hand, have so far been utilized in two ways (Caserta & Voß 2010):

• •

as general-purpose optimization methods ready to apply to any problem without any modification, and

• •

as higher-order methods which guide how problem-specific heuristics are applied to instances belonging to a particular problem class.

Over time, it has been noticed that the latter approach yields better results (Caserta & Voß 2010). However, once a researcher evaluates such a method, they assess the combined performance of the metaheuristic and its subordinate heuristic(s). Although that gives insight into performance of the method as a whole, it does not provide answers to the following questions:

• •

Is it possible that the performance score has been achieved mostly or solely by the heuristic(s)?

• •

How much does the guiding logic of the metaheuristic contribute to total performance?

The answers to these questions are important because if it is the case that performance comes mostly or solely from heuristics, then it would be wrong to attribute the score to the metaheuristic and claim that a new solver for the specific class of problems has been found. The goal of this paper is to present a sound methodological framework to answer said questions. To our best knowledge, this is the first attempt to formulate such a technique.

The rest of the paper is organized as follows. The proposed methodology is described in Section 2. In Section 3, we present an example of its application to Simulated Annealing, Boolean Satisfiability Problem and the Flip heuristic. Finally, we discuss it and draw our conclusions in Section 4.

2 Proposed Methodology

Let $\mathcal{M}$ be the metaheuristic being examined, and let $\mathcal{H}$ denote a single heuristic or a group of heuristics $\mathcal{H}=\{\mathcal{H}_{1},\mathcal{H}_{2},\ldots,\mathcal{H}_{m}\}$ intended to be executed one after another, known to work well on the problem class of interest. The performance metric of the combined method $\mathcal{M}[\mathcal{H}]$ in which $\mathcal{M}$ guides the application of $\mathcal{H}$ can be modeled as a random variable $Y$ whose distribution is given by:

[TABLE]

where $\Pi$ and $S$ are random variables representing the instance of the problem class to which $\mathcal{M}[\mathcal{H}]$ is applied and the seed for the random number generator, whereas $\theta_{\mathcal{M}}$ and $\theta_{\mathcal{H}}$ denote the parameters of $\mathcal{M}$ and $\mathcal{H}$, respectively. For now, we assume that $Y$ is a univariate variable, i.e. that the metric is a single value (the objective function to optimize, execution time, etc.). Its distribution is not known in advance and researchers estimate it by first tuning $\theta_{\mathcal{M}}$ and $\theta_{\mathcal{H}}$ and then evaluating the method on a number of problem instances $\pi_{1},\pi_{2},\ldots,\pi_{l}$, repeating evaluation several times for different choices of the seeds for the random number generator.

As said in introduction, the metric $Y$ measured in this manner represents an estimate of the performance of $\mathcal{M}[\mathcal{H}]$. In order to assess how good $\mathcal{M}$ is at guiding $\mathcal{H}$, we can introduce an additional variable $M$ to Equation 1 which now becomes:

[TABLE]

assuming a more general form for the performance of a metaheuristic ($M$) guiding $\mathcal{H}$ for the problem of a given class. The variable $M$ will denote the metaheuristic component and will be understood to have two levels: $\mathcal{M}$ and $\emptyset$, where the latter denotes what we will call a naive or placebo metaheuristic henceforth. It is a metaheuristic which is based on no purposeful logic and has no components other than random decisions. It is such that $\emptyset[\mathcal{H}]$ acts as an algorithm where $\mathcal{H}$ is guided randomly, as if no metaheuristic has been used to guide it. Then, to answer the question:

• •

How good is $\mathcal{M}$ at guiding $\mathcal{H}$?

we should estimate

[TABLE]

and compare it to:

[TABLE]

The difference reveals the effect of changing the metaheuristic from naive random search, which has no guiding logic, to $\mathcal{M}$. If the effect is negligible, then it indicates that using $\mathcal{M}$, which may have sophisticated and complicated logic, to guide how $\mathcal{H}$ is applied, is the same as using a naive metaheuristic with no logic to guide $\mathcal{H}$. In fact, that would mean that any score $\mathcal{M}[\mathcal{H}]$ has achieved comes from using $\mathcal{H}$ and has nothing or little to do with $\mathcal{M}$. After all, if the logic of $\mathcal{M}$ guides $\mathcal{H}$ similarly or identically to random search, then we cannot justify use of $\mathcal{M}$ in that particular setting.

This method is similar to the one used in a typical scenario where there are two factors, $A_{1}$ and $A_{2}$, and a researcher wants to estimate the linear effect of $A_{1}$ on a yield variable $Y$ when $A_{2}$ is fixed to a certain value. The way to do so is to define the low and high levels of $A_{1}$ and then estimate how $Y$ changes when $A_{1}$ is increased from its low to high level. The effect that we are estimating is called the simple effect of $A_{1}$ at the chosen level of $A_{2}$. This is precisely what we are trying to do in our case. We want to estimate how the performance metric changes when $M$, the metaheuristic component, is changed from its low level with no logic ($\emptyset$), to its high level, the metaheuristic $\mathcal{M}$ being examined, with the heuristic component fixed to $\mathcal{H}$. The method that we propose in this Section achieves just that. Another example analogous to our case is from pharmacological studies. When a new medicine is tested, one group of patients, called the control group, is given placebo, while the other is given the drug. If the effect of the medicine is significantly better then that of placebo, the drug is deemed effective. Otherwise, there is no justification to produce and use the medicine as it is less effective than a simple placebo. If we use terminology from that example in our study, we will say that the $\emptyset$ acts as placebo and $\mathcal{M}$ takes the role of medicine.

We have exposed the core of our methodology and the rationale behind it, but there are still several issues that we must address:

Is $\emptyset$ the same for every $\mathcal{M}$ and why cannot we simply apply $\mathcal{H}$ without plugging it in a metaheuristic? 2. 2.

Should the parameters be tuned to yield the maximal performance prior to evaluation of $\mathcal{M}[\mathcal{H}]$ and $\emptyset[\mathcal{H}]$ or drawn randomly from a predefined space of allowed values? 3. 3.

How to compare $P(Y|M=\mathcal{M},\Pi,S,\theta_{\mathcal{M}},\theta_{\mathcal{H}})$ to $P(Y|M=\mathcal{\emptyset},\Pi,S,\theta_{\emptyset},\theta_{\mathcal{H}})$?

We answer all those questions in continuation of this section.

2.1 The Naive Metaheuristic

The main idea of our method is to see if guiding $\mathcal{H}$ with no logic is the same as guiding it with the logic of $\mathcal{M}$, the metaheuristic being examined. The rationale behind this is that each metaheuristic is a specific set of rules, and that if using those rules gives the same results as not using any rules at all, then the observed performance is achieved by heuristics alone and the logic of $\mathcal{M}$ is not effective nor efficient. We have referred to guiding heuristics with no logic as the naive or placebo metaheuristic, $\emptyset$. The reason why, in general, $\emptyset$ has to be a naive metaheuristic, and not just a mere application of $\mathcal{H}$ is that $\emptyset[\mathcal{H}]$ has to invest the same computational effort as $\mathcal{M}[\mathcal{H}]$ in order for comparison of the corresponding distributions to be fair. This means that if $\mathcal{M}$ is a population metaheuristic (such as, e.g., Genetic Algorithm), $\emptyset$ must be too. Similarly, if $M$ is a single solution metaheuristic (such as, e.g., Simulated Annealing), so must be $\emptyset$. Moreover, in the former case, if the population in $\mathcal{M}$ consists of $N$ individuals, the same must hold for $\emptyset$.

In general, we do not need to state $\emptyset$ explicitly. We can derive $\emptyset[\mathcal{H}]$ from $\mathcal{M}[\mathcal{H}]$ by removing all $\mathcal{M}$’s unique algorithmic components and leaving only naive, random operations. An example in Section 3 will clarify this step.

2.2 Choice of Parameters

The choice of parameters is crucial to performance of a (meta)heuristic. If tuned appropriately, they can greatly improve performance. If not set to appropriate values, they can deteriorate the algorithm. The question that naturally arises in our case is whether the parameters should be tuned prior to evaluation or treated as random variables and randomly set before each run of $\mathcal{M}[\mathcal{H}]$ and $\emptyset[\mathcal{H}]$ during their evaluation. Both alternatives are viable, but are related to essentially different goals. If we opt for randomly setting the parameters, we would be aiming to assess the intrinsic guiding capability of $\mathcal{M}$ that does not depend on the choice of parameters and is present in all its applications. However, does such capability exist? Different parameter settings can lead to diametrically opposite results. Besides, before $\mathcal{M}[\mathcal{H}]$ is applied to real problems in practice, it is always tuned. Practitioners and researchers are interested in the best performance $\mathcal{M}[\mathcal{H}]$ can give for a class of problems, not any performance for random parameter settings. Hence, we argue for tuning the parameters of $\mathcal{M}[\mathcal{H}]$ prior to its evaluation. We can use some of them as the parameters of $\emptyset[\mathcal{H}]$ (for example, the number of individuals to ensure the populations are of the same sizes in $\mathcal{M}$ and $\emptyset$), and then tune the parameters of $\mathcal{H}$, if any.

2.3 Comparison of Distributions

In literature, the most common way to compare two metaheuristics is to compare their expected values of $Y$, the chosen performance metric, approximated by means of measurements of $Y$ on the selected problems instances for different, but random choices of the seed for the random number generator. However, we argue against using means to compare distributions of $Y$. What must be understood is that mean, even when accompanied by standard deviation, may not be representative of the distribution (Gunawardena 2014). Therefore, difference in means may not be informative and inference based on it may be invalid. Another, unfortunately common practice that we argue against is using $p$ values as definite proofs to accept or reject tested hypotheses. One reason is that significance at the desired level can always be achieved by using sufficiently large samples (Demidenko 2016). In our case, by evaluating algorithms on a large number of problem instances and repeating the process for a lot of times, we can make $p$ values as small as desired. The other reason is, as Fraser & Reid (2016) explained, that ”[$p$ value] can guide the judgments about scientific conclusions, but cannot replace them.”.

Knowing this, we ask what is the appropriate way to compare the distribution of the performance metric for $\mathcal{M}[\mathcal{H}]$, $Y_{\mathcal{M}[\mathcal{H}]}$, with that for $\emptyset[\mathcal{H}]$, $Y_{\emptyset[\mathcal{H}]}$? Let us assume, without loss of generality, that lower scores of the metric signify superior performance. If $\mathcal{M}[\mathcal{H}]$ works better than $\emptyset[\mathcal{H}]$, we should expect the distribution of $Y_{\mathcal{M}[\mathcal{H}]}$ to be located to the left of $Y_{\emptyset[\mathcal{H}]}$. A measure of how far the former is to the left of the latter is

[TABLE]

the probability that a score of a run of $\mathcal{M}[\mathcal{H}]$ is lower, i.e. better than that of a score of a run of $\emptyset[\mathcal{H}]$ for some randomly selected problem instance. However, we should not limit ourselves to testing only if $Y_{\mathcal{M}[\mathcal{H}]}$ is located to the left of $Y_{\emptyset[\mathcal{H}]}$. For example, if the scores of $\mathcal{M}[\mathcal{H}]$ lied in the range $(1.1\times 10^{-4},1.2\times 10^{-4})$ and those of $\emptyset[\mathcal{H}]$ in $(1.3\times 10^{-4},1.4\times^{-4})$, but score differences lower than $10^{-3}$ are practically negligible, even though $P(Y_{\mathcal{M}[\mathcal{H}]}<Y_{\emptyset[\mathcal{H}]})$ would be equal to $1$ and indicate complete superiority of $Y_{\mathcal{M}[\mathcal{H}}$ over $Y_{\emptyset[\mathcal{H}]}$, which would be true from a purely statistical point of view, but false from the standpoint of practical importance. Therefore, we first need to set some threshold $\delta\geq 0$ to define a minimal difference between two performance scores required for them to be considered practically different. Therefore, instead of estimating $P(Y_{\mathcal{M}[\mathcal{H}]}<Y_{\emptyset[\mathcal{H}]})$, we should focus on:

[TABLE]

The probability quantifies practical benefit, with respect to $\delta$, of guiding $\mathcal{H}$ with $\mathcal{M}$, hence the name $b$. The converse probability

[TABLE]

represents the risk of using $\mathcal{M}[\mathcal{H}]$ instead of $\emptyset[\mathcal{H}]$, that is the probability that $\mathcal{M}$ guides $\mathcal{H}$ practically worse than $\emptyset$. What remains is the probability that $\mathcal{M}[\mathcal{H}]$ and $\emptyset[\mathcal{H}]$ are practically equivalent:

[TABLE]

Those quantities, which we will call BER values (benefit, equivalence, risk) from now onwards, express the size of the effect of using $\mathcal{M}$ instead of $\emptyset$ on the probability scale, simultaneously taking into account chosen definition of practical meaningfulness. The BER values are related to ROC curves Gonçalves et al. (2014). More specifically, when $\delta=0$, the $b$ is the area under the ROC curve (AUROC) associated with $Y_{\mathcal{M}[\mathcal{H}]}$ and $Y_{\emptyset[\mathcal{H}]}$, whereas $r+e$ is equal to the area above the curve (Demidenko 2016). This is not a new idea for testing for difference between two distributions. We refer interested readers to (Wolfe & Hogg 1971, Zhou 2008, Newcombe 2006a, b, Demidenko 2016) for more details about ROC curves, computational techniques for estimating AUROC, and application of the method to discriminate between distributions. What is new in our approach is $\delta$, the threshold of practical significance, which should be set in advance according to the theory and empirical knowledge of the optimization-problem class for which $\mathcal{M}[\mathcal{H}]$ is being developed, and making distinction between $r$ and $e$ values - as opposed to Demidenko (2016) who does not distinguish between them.

Finally, we have to address calculation of $b$, $r$, and $e$ and their interpretation. Let us assume that we have run $\mathcal{M}[\mathcal{H}]$ and $\emptyset[\mathcal{H}]$ on problem instances $\pi_{1},\pi_{2},\ldots,\pi_{l}$, repeating evaluation on each instance $n$ times using seeds $s_{ij}$, $i=1,2,\ldots,l$, $j=1,2,\ldots,n$. Let $\mathbf{Y}^{*}_{\mathcal{M}}$ and $\mathbf{Y}^{*}_{\emptyset}$ denote $l\times n$ matrices where the results of $\mathcal{M}[\mathcal{H}]$ and $\emptyset[\mathcal{H}]$ are stored. The obvious way to calculate empirical $b$ value, denoted as $b^{*}$, is to compare the corresponding entries in the result matrices:

[TABLE]

where $I_{\varphi}$ is the indicator function that takes the value $1$ when its underlying condition $\varphi$ evaluates to $\top$, and [math] otherwise.

The explanation of Equation 9 is as follows. The result of the sums in Equation 9 is equal to the number of times that $\mathcal{M}[\mathcal{H}]$ produced better solutions than $\emptyset[\mathcal{H}]$ for the problems $\pi_{1},\pi_{2},\ldots,\pi_{l}$. The denominator $l\times n^{2}$ is the total number of comparisons. Therefore, their ratio is an estimate of the probability that for a random instance from the problem class to which $\pi_{1},\pi_{2},\ldots,\pi_{l}$ belong, $\mathcal{M}[\mathcal{H}]$ will produce a better solution than $\emptyset[\mathcal{H}]$. Better in this context means ”lower for at least $\delta$”.

The empirical equivalence and risk are calculated analogously:

[TABLE]

[TABLE]

In general, if $\mathcal{M}$ is a good choice for guiding $\mathcal{H}$, we should expect that $b^{*}>1/2$ and $r^{*}\approx 0$. If $e^{*}\approx 1$, then no meaningful difference between $\mathcal{M}[\mathcal{H}]$ and $\emptyset[\mathcal{H}]$ has been found, suggesting that the score of $\mathcal{M}[\mathcal{H}]$ is achieved by $\mathcal{H}$. The greater the value of $r^{*}$, the stronger the evidence that $\mathcal{M}$ guides $\mathcal{H}$ in a way that it deteriorates the effect of the heuristic component. The closer $b^{*}$ is to $1$, the stronger the evidence in favor of $\mathcal{M}$ being able to efficiently guide $\mathcal{H}$.

Finally, we must stress out that BER values, just as $p$ value, cannot replace a researcher’s own reasoning. Just as we should not base conclusions solely on $p$ values being lower or greater than the usual significance thresholds of $0.01$ and $0.05$, we should not regard empirical benefit, risk, and equivalence as a definite answer to the question concerning the examined metaheuristic’s efficiency in guiding its heuristic(s). After all, the nature of statistical research is such that only through replications of experiments can a certain hypothesis be accepted or rejected. So, researchers should always plot $\mathbf{Y}^{*}_{\mathcal{M}}$ and $\mathbf{Y}^{*}_{\emptyset}$ one against another to visually inspect the empirical distributions. Moreover, similar plots should be made for each problem instance. Only when all that is taken into account, should researchers formulate their conclusions.

2.4 Assumptions

We will conclude this section by briefly stating the assumptions of the methodology which we proposed:

1. A1

The heuristic component $\mathcal{H}$ is known to work well. 2. A2

The performance metric $Y$ is univariate (its value is a single score, not a tuple). 3. A3

The metric $Y$ is measured for each run.

We can see that they are fairly general and easy to meet in practice. In Section 4, we discuss the cases of their violation.

3 Experimental Example

In this Section, we will describe how we applied the method presented in Section 2 to a variant of Boolean Satisfiability Problem, 3-SAT, Simulated Annealing (SA), and the SAT heuristic known as Flip. The problem and the algorithms are presented in Sections 3.1-3.6. We describe benchmarks in Section 3.7, tuning in Section 3.8, and the results of comparing SA[Flip] to $\emptyset$[Flip] in Section 3.9.

The repository with the code and data can be downloaded from https://osf.io/f2m9w/. .

3.1 Boolean Satisfiability Problem

Boolean Satisfiability Problem, shorthand SAT, is an NP-complete problem (Cook 1971) formulated as follows:

Definition 1

Given a Boolean formula $F$ with $n$ propositional letters $p_{1},p_{2},\ldots,p_{n}$, find their valuation under which $F$ evaluates to $\top$.

Any Boolean formula can be converted to conjunctive normal form (CNF), i.e. a conjunction of clauses that are themselves disjunctions of literals (propositional letters or their negations):

[TABLE]

where $\pm p_{i,j}\in\{p_{l},\neg p_{l}\}$ for some $l\in\{1,2,\ldots,n\}$. If all $k_{i}=k$, then we say that $F$ is in its $k$-CNF and refer to the SAT problem as $k$-SAT. Since $3$-SAT is also NP-complete and all Boolean formulae can be converted to $3$-CNF (Cook 1971), we will focus on the case where $k=3$.

3.2 Solution Representation and Objective Function

A solution to a ($3$-)SAT instance is a valuation of its propositional letters $p_{1},p_{2},\ldots,p_{n}$, i.e. a mapping from $\{p_{1},p_{2},\ldots,p_{n}\}$ to $\{\bot,\top\}$. Encoding $\bot$ as [math] and $\top$ as $1$, we can represent solutions as integer arrays of zeros and ones. The goal is to find such a solution that all the clauses in the formula are satisfied. We can formulate the objective function as the percentage of satisfied clauses and aim to maximize it, or the percentage of unsatisfied clauses and try to minimize it. The two objective functions are equivalent and both return the values between [math] and $1$. We opt for the minimization alternative in this paper (i.e. we will minimize the ratio of the number of unsatisfied clauses to the total number of clauses) and use it as the performance metric $Y$. Obviously, if $Y=0$, the optimal solution has been found and satisfiability of the formula in question has been proven. The reason why we use percentages rather than numbers of unsatisfied clauses is that we want the performance metric to be on the same scale for all formulae, no matter how much clauses they consist of.

3.3 Flip Heuristic

Let $F$ be a $3$-CNF formula. The Flip heuristic receives a possible solution $\mathbf{x}=[x_{1},x_{2},\ldots,x_{n}]$ ($x_{i}\in\{0,1\}$ for $i=1,2,\ldots,n$) and iteratively flips one its element at a time if it improves the objective function until no further improvement is possible. The heuristic is presented in Algorithm 1 (Marchiori & Rossi 1999).

3.4 Simulated Annealing

Simulated Annealing (SA) is a well-known and widely used metaheuristic whose history dates back to 1980s when Kirkpatrick et al. (1983) and Černỳ (1985) published first papers on the algorithm. Simulated Annealing is inspired by the process of physical annealing with solids, ”in which a crystalline solid is heated and then allowed to cool very slowly until it achieves its most regular possible crystal lattice configuration (i.e., its minimum lattice energy state), and thus is free of crystal defects.” (Nikolaev & Jacobson 2010). The algorithm starts with an initial solution and processes it iteratively. Each iteration consists of several steps, and at each step, the algorithm compares the current solution to one if its neighbors. The current solution is always replaced with the better neighbor. If the neighbor is worse, replacements occur with a probability which depends on the current temperature. The algorithm receives the initial temperature and the cooling schedule at the beginning and decreases temperature at each iteration according to the schedule. The pseudo-code of the Simulated Annealing is outlined in Algorithm 2 (Nikolaev & Jacobson 2010).

3.5 The SA[Flip] Algorithm for the SAT problem

The combination of Simulated Annealing and the Flip heuristic for SAT, named SA[Flip], is presented in Algorithm 3. In it, the heuristic specific to SAT is applied to the initial solution at the beginning of the algorithm, and once to each neighbor proposed to replace the current solution. Even though there may be other ways to combine the two algorithms, the goal of our study is not to find the best combination of them all, but to show how we can assess if the overall result of the combination being examined is due to the heuristic alone. The same procedure can be carried out for any metaheuristic and the heuristic(s) it guides.

We used geometric cooling schedule, in which $T_{k+1}=\alpha T_{k}$ for some constant $\alpha\in(0,1)$. As the stopping criterion we used the following compound condition:

• •

The objective function (ratio of the number of unsatisfied clauses to the total number of clauses) of the current solution is equal to [math] or

• •

$k$, the number of iterations performed, is equal to $MNI$, the maximal number of iterations allowed, specified as a SA[Flip]’s parameter.

We checked for the stopping condition at each iteration as well as after each step.

Also, we kept track of the best solution encountered during execution of the algorithm and output it when SA[Flip] stops. We decided to do so because it may happen that the algorithm finds the optimal solution, but replaces it with a neighbor that is worse than it.

Even though definition of neighborhoods can be treated as an additional parameter to calibrate, we chose not do so, but to adopt one neighborhood definition in advance in order to reduce the number of parameters and simplify demonstration of our methodology. Of course, we advise researchers to experimentally determine the best definition of a neighborhood, as in (Simić 2017). The one that we adopted and used throughout the experiment is as follows:

Definition 2

Two solutions to the same instance of $3$-SAT problem are neighbors to each other if and only if their Hamming distance is equal to $1$.

This means that a neighbor of a solution differs from it in valuation of a single propositional letter.

3.6 Derivation of $\emptyset$[Flip]

As said in Section 2, we need to compare SA[Flip] to $\emptyset$[Flip] in order to estimate how good SA is at guiding Flip. We do not need to explicitly state $\emptyset$ as an actual algorithm. It is sufficient to remove all SA’s components from SA[Flip] and leave only naive operations at the metaheuristic level: random generation of the initial solution, random generation of neighbors, and their random acceptance. The parameters inherited from SA[Flip] are $MNI$ and $M$ and they should be set to the same values as for SA[Flip] in order to ensure that $\emptyset$[Flip] can invest the same computational effort as SA[Flip].

We present $\emptyset$[Flip] in Algorithm 4.

3.7 Benchmarks

Even though $3$-SAT constitutes a class of problems of its own, we did not aim to cover all the possible subclasses of $3$-SAT problems. Instead, we focused on those $3$-SAT instances which are in the so called phase transition. Those are the formulae with approximately $4.24n$ clauses (Gent & Walsh 1994), where $n$ is the number of propositional letters that appear in them. Such instances are computationally hardest to solve and the probability of them being satisfiable is approximately equal to the probability that they are not.

We also limited $n$ since it is impossible to conduct an experiment involving all possible numbers of propositional letters and our computational resources were limited. We chose the range $50\leq n\leq n$ because the corresponding solution spaces are sufficiently large but not too much for our testing machine. For $n=50,75,100,125$, we downloaded $100$ corresponding instances from SATLIB (http://www.cs.ubc.ca/~hoos/SATLIB/benchm.html) (Hoos & Stützle 2000). They are all satisfiable and in the phase transition. We split the formulae into training and test sets. The former contained $100$ formulae, $20$ for each $n=50,75,100,125$, while the latter included the rest.

3.8 Tuning the Parameters

We tuned the parameters following the methodology of Simić (2017) as it rigorously employs statistical techniques from Design of Experiments (Montgomery 2000).

First, we screened the parameters $T_{0}$, $\alpha$, $M$, and $MNI$ to identify the influential ones. To do so, we defined their low, medium, and high levels (see Table 1). We used Box-Behnken design for four three-level factors (Oehlert 2000, Box & Behnken 1960) and evaluated SA[Flip] for thirty times on each formula in the training set, blocking the design for seeds. It turned out that $M$ and $MNI$ had substantial main effects and that there were second-order interactions between $T_{0}$ and $\alpha$, on one hand, and $M$ and $MNI$ on the other. Hence, we had to calibrate all the four parameters. The found effects are presented in Figure 1.

Then, we calibrated the parameters iteratively, conducting Response Surface Methodology #, evaluating SA[Flip] for thirty times on each benchmark, but without blocking the design for seeds. The design that we used in this phase was $2^{4-1}$ fractional factorial. The reason why we used such a simple design is that the response (average performance) can be approximated with a linear model if the portion of the search space is sufficiently small. To ensure that, we used small but effective half-distances for the parameters. We present them in Table 2. The starting configuration was: $T_{0}=50$, $\alpha=0.9$, $M=20$ and $MNI=50$, because screening indicated that it might give very good results. We stopped the procedure once the values of $M$ and $MNI$ were such that the maximal number of applications of Flip exceeded $5000$. The found settings are: $T_{0}=51.71$, $\alpha=0.92$, $M=50$ and $MNI=103$.

3.9 Results and Their Interpretation

We evaluated SA[Flip] and $\emptyset$[Flip] on each testing benchmarks thirty times. We made sure that the algorithms used the same sets of seeds for each formula to allow for fair comparison.

In general, both methods achieved very good results, as can be seen in Tables 3-6. Their performance scores, $Y$, deteriorate as $n$ increases, which is also the case with their success rates - the percentages of successful runs, i.e. the runs where the output $Y=0$.

We calculated BER values for each $n=50,75,100,125$ as well as for the whole test set. We used three different values for $\delta$ in our analysis: [math], $0.01$, and $0.02$. The results are presented in Tables 7-9 and depicted in Figure 2.

Overall, the $e$ value turned out to dominate other two by large margins for all choices of $\delta$ and each $n=50$, $75$, $100$, and $125$. This implies that SA[Flip] is effectively the same as $\emptyset$[Flip], i.e. that SA guides Flip as effectively as the corresponding naive metaheuristic. We can also observe that $e^{*}$ drops whereas $b^{*}$ and $r^{*}$ increase with $n$ for $\delta=0$. It is probable that such a trend continues for $n>125$ and SA[Flip] and $\emptyset$[Flip] become effectively distinct at some point. Therefore, future research could focus on investigating this hypothesis.

Plots of distributions of $Y$ for SA[Flip] and $\emptyset$[Flip] are presented in Figure 3. By visual inspection, we conclude that the distributions are almost indistinguishable, which confirmes what $e^{*}$ has indicated: that SA is not than $\emptyset$ at guiding Flip (for this set of problems). In turn, that means that the observed performance of SA[Flip] is most probably due to Flip alone. Had we not tested SA[Flip] in this manner, we would not have discovered that efficiency of SA[Flip] came from the heuristic component. We were able to find it out only because we compared SA[Flip] to $\emptyset$[Flip], which stresses out the importance of the methodology proposed in this paper and its usefulness in research in this field.

There is one more issue that needs to be discussed. Another explanation for observed equivalence of SA[Flip] and $\emptyset$[Flip] could be that the benchmarks we used are too easy so both algorithms performed really well and no difference between them was possible to found in the first place. This is known as ceiling effect (Bartz-Beielstein & Preuß 2014). However, we do not think that the effect occurred in this experiment because we deliberately used the benchmarks that are in phase transition and hence, the most challenging and difficult to solve. In addition to this, the numbers of propositional letters were not low and there are studies which evaluated solvers on the same groups of benchmarks but reported worse results, e.g. (Djenouri et al. 2016). Still, researchers who decide to follow our methodology to estimate their metaheuristics should take caution and make sure that their benchmarks do not cause ceiling (or floor) effects.

Finally, our conclusion is as follows. Success of our SA[Flip] for the class of random $3$-SAT formulae in phase transition with $125$ propositional letters at most comes from the Flip heuristic. For this class of Boolean formulae, guiding Flip with Simulated Annealing has the same effect as using a naive random metaheuristic with essentially no logic to guide application of Flip. The more similar SAT problems are to those used in our study, the higher is the probability that the same effects will be detected.

4 Discussion and Conclusions

In this paper, we have proposed a methodology to empirically estimate how efficient a metaheuristic algorithm $\mathcal{M}$ is at guiding specific heuristic(s) $\mathcal{H}$. The proposed technique was applied to Simulated Annealing (SA), Boolean Satisfiability Problem and the Flip heuristic. The experiment revealed that the performance score of the combination of SA and Flip was due to the heuristic, which is a result that we would not be able to obtain without our methodology.

The methodology itself is mathematically well-founded, intuitive, easy to apply and relies on practical significance rather than solely on statistics. It directly compares empirical distributions of the chosen performance metrics, not just sample means, which provides a better insight into the metaheuristic in question. By comparing $\mathcal{M}[\mathcal{H}]$ to $\emptyset[\mathcal{H}]$, the technique allows us to estimate the effect of using metaheuristic $\mathcal{M}$ to guide $\mathcal{H}$. The introduced BER values quantify that effect on the probability scale and, accompanied by visual comparison, reveal whether it is justified to guide application of $\mathcal{H}$ with $\mathcal{M}$. Without investigating if performance of $\mathcal{M}[\mathcal{H}]$ comes mostly or entirely from $\mathcal{H}$, we can easily draw wrong conclusions and claim that we have discovered novel solvers, when, if fact, we have done nothing more than wrapping up an efficient heuristic solver with a metaheuristic whose contribution is negligible. If we propose $\mathcal{M}[\mathcal{H}]$ as a new solver, we must prove that there is something that makes it worth to guide $\mathcal{H}$ with $\mathcal{M}$, and the technique studied and demonstrated in this paper offers a way to do precisely that.

Comparison to $\emptyset[\mathcal{H}]$ is the core of our approach. We defined $\emptyset$ to be a naive, placebo metaheuristic, which performs only naive operations in guiding $\mathcal{H}$. We argued that completely random decisions constitute naive moves and that such $\emptyset$ represents the ”low” level of $M$ in Equation 2, equivalent to no guiding logic. Can there be naive operations other than completely random decisions? Can greedy moves be thought of as naive? Those would be operations that always select the best solution from a group of candidates and discard the rest. Although they may seem naive, they do follow some logic, no matter how simple it is. Therefore, $\emptyset[\mathcal{H}]$ that would include greedy moves would actually be a plain greedy algorithm. If $\mathcal{M}[\mathcal{H}]$ fails to beat it, we can say that there are no reasons to use $\mathcal{M}[\mathcal{H}]$ when a simple solver achieves the same or better results, but we would not be able to test if performance of $\mathcal{M}[\mathcal{H}]$ is achieved by $\mathcal{H}$ or the logic of $\mathcal{M}$ contributes to it significantly. Therefore, we argue for $\theta[\mathcal{H}]$ to contain only random operations.

One may also ask why we do not compare $\mathcal{M}[\mathcal{H}]$ to $\mathcal{M}[\emptyset]$, where $\emptyset$ would denote use of no heuristic at all? The reason is that difference between $\mathcal{M}[\mathcal{H}]$ and $\mathcal{M}[\emptyset]$ can reveal if $\mathcal{H}$ contributes anything to the performance of the whole method, not if $\mathcal{M}$ is able to guide it efficiently.

The proposed approach is not without limitations, though. First of all, it requires evaluation of an additional algorithm ($\emptyset[\mathcal{H}]$) which is derived from the metaheuristic being examined. Even though this prolongs research, it also provides information which we would not get otherwise, as demonstrated in the example in Section 3, and without which scientific conclusions could be flawed. Therefore, we find that taking more time to complete this step pays off. Another limitation is that it enables us to reason about our metaheuristic only with respect to a chosen class of optimization problems. However, this limitation is not unique to this methodology and is inherent to all techniques for analyzing numerical experiments involving stochastic optimization algorithms.

Also, we must ask if the hybrid algorithm $\mathcal{M}[\mathcal{H}]$ is the only (sensible) combination of $\mathcal{M}$ and $\mathcal{H}$ because if it is the case, then we can be sure that our method sheds light onto general ability of $\mathcal{M}$ to guide $\mathcal{H}$. However, there may be more than one way to guide $\mathcal{H}$ with $\mathcal{M}$. For instance, had we used Genetic Algorithm (GA) (Holland 1992) instead of SA in the example in Section 3, we could have applied Flip after mutation (as Marchiori & Rossi (1999)), but we could have also done it before mutation, immediately after performing crossovers. In such cases, rather than estimating general ability of $\mathcal{M}$ to guide $\mathcal{H}$, we are assessing efficiency of the specific strategy based on $\mathcal{M}$ for guiding $\mathcal{H}$. If its effect is approximately equal to that of $\emptyset$, which can be tested with our method, then we can determine if guiding $\mathcal{H}$ in that particular way is justified. Moreover, even though it is possible to plug a heuristic into a metaheuristic between any two operations, is it really sensible to arbitrary intertwine their logics? In each metaheuristic there is a point where the quality of a solution, i.e. the value of the objective function, is computed. It is the only step in execution of metaheuristic methods where they need to evaluate a problem-specific function. All the other operations that they perform are based on specific optimization ideas or some metaphors and constitute a logical unity. In our opinion, the moment just before evaluating a solution is a good time for applying a heuristic because that keeps problem-specific operations (evaluation of the objective function and application of heuristic) at one place, allowing the logic of the metaheuristic to execute without interruptions and as originally designed. This way, comparing $\mathcal{M}[\mathcal{H}]$ to $\emptyset[\mathcal{H}]$ is as close to revealing the general effect of using $\mathcal{M}$ to guide $\mathcal{H}$ as it gets.

We also need to discuss the assumptions of our methodology as well as the effects of violating them. One of the assumptions is that $\mathcal{H}$ is an efficient heuristic. If $\mathcal{H}$ is new and its efficiency has not been confirmed, then the heuristic must be tested prior to application of the proposed methodology. Another assumption is that performance metric $Y$ is univariate, i.e. a single value, not a tuple of values. If several metrics are of interest, we can compare $\mathcal{M}[\mathcal{H}]$ to $\emptyset[\mathcal{H}]$ once for each of them and then analyze results per metric. Then, what if we want to use a metric calculated as an aggregated value of the results of several runs on each problem instance? Let us suppose that we have stored the results in $l\times n$ matrices $\mathbf{Y}_{\mathcal{M}}^{*}$ and $\mathbf{Y}_{\emptyset}^{*}$. The actual metric scores that we are interested in are then calculated as $Z_{\mathcal{M}}^{*}[i]=f(\mathbf{Y}^{*}_{\mathcal{M}}[i,1],\mathbf{Y}^{*}_{\mathcal{M}}[i,2],\ldots,\mathbf{Y}^{*}_{\mathcal{M}}[i,n])$ and $Z_{\emptyset}^{*}[i]=f(\mathbf{Y}^{*}_{\emptyset}[i,1],\mathbf{Y}^{*}_{\emptyset}[i,2],\ldots,\mathbf{Y}^{*}_{\emptyset}[i,n])$ ($i=1,2,\ldots,l$), where $f$ is the aggregating function. The methodology can still be applied, but the formulae for empirical BER values would need to be modified and the results could not be interpreted in quite the same way. Instead of comparing $P(Y|M=\mathcal{M},\Pi,S,\theta_{\mathcal{M}},\theta_{\mathcal{H}})$ to $P(Y|M=\mathcal{\emptyset},\Pi,S,\theta_{\emptyset},\theta_{\mathcal{H}})$, we would essentially be comparing $P(Z|M=\mathcal{M},\Pi,\theta_{\mathcal{M}},\theta_{\mathcal{H}})$ to $P(Z|M=\mathcal{\emptyset},\Pi,\theta_{\emptyset},\theta_{\mathcal{H}})$ and the value $b=P(Z_{\mathcal{M}}<Z_{\emptyset}-\delta)$ would answer the following question:

• •

What is the probability that, for a randomly chosen instance from the problem class of interest, $\mathcal{M}[\mathcal{H}]$’s score will be practically better than that of $\emptyset[\mathcal{H}]$ when aggregated over several runs?

This is different from the meaning of the $b$ value as originally defined in Section 2 for the non-aggregated case. The corresponding formula for $b^{*}$ would then be:

[TABLE]

with analogous modifications being in place for $e^{*}$ and $r^{*}$. Those differences are simple, but subtle, so we need to point them out.

Finally, the BER values that we define and propose to quantify the degree to which two distributions are not just statistically, but practically different, can be used to compare any two stochastic algorithms, not just $\mathcal{M}[\mathcal{H}]$ and $\emptyset[\mathcal{H}]$. Moreover, since numerical and practical significance are confounded in BER values through $\delta$, we find them suitable to detect important effects not just in the field of metaheuristics, but in science in general.

We hope that other researchers will see merit in our idea, adopt it in their own studies and improve it further to the benefit of the whole research community.

Possible directions of future research are:

• •

Developing a methodology that would simultaneously test both the efficiency of $\mathcal{H}$ and $\mathcal{M}$’s ability to guide it;

• •

Formulating a technique capable of estimating the general ability of $\mathcal{M}$ to guide any heuristic for the problem at hand, not just the selected $\mathcal{H}$.