Exploiting Problem Structure in Combinatorial Landscapes: A Case Study on Pure Mathematics Application

TL;DR

This paper introduces an AI-based method that leverages problem structure in combinatorial landscapes to efficiently solve a pure mathematics problem involving narrow admissible tuples, demonstrating significant search space reduction and solution effectiveness.

Contribution

It presents a novel approach that exploits local search structures to improve combinatorial optimization in pure mathematics applications.

Findings

Efficiently finds best known solutions

Reduces search space significantly

Effectively escapes local minima

Abstract

In this paper, we present a method using AI techniques to solve a case of pure mathematics applications for finding narrow admissible tuples. The original problem is formulated into a combinatorial optimization problem. In particular, we show how to exploit the local search structure to formulate the problem landscape for dramatic reductions in search space and for non-trivial elimination in search barriers, and then to realize intelligent search strategies for effectively escaping from local minima. Experimental results demonstrate that the proposed method is able to efficiently find best known solutions. This research sheds light on exploiting the local problem structure for an efficient search in combinatorial landscapes as an application of AI to a new problem domain.

Equations

$$H(k)\le k\log k+k\log \log k-(1+\log 2)k+o(k).$$

$$H(k)\le k\log k+k+o(1).$$

Full text

1 Introduction

AI techniques have shown their advantages on solving different combinatorial optimization problems, such as satisfiability Sutton et al. (2010); Dubois and Dequen (2001); Bjorner and Narodytska (2015); Ansótegui et al. (2015), traveling salesman problem Zhang (2004), graph coloring Culberson and Gent (2001), job shop scheduling Watson et al. (2003), and automated planning Bonet and Geffner (2001).

These problems can be generalized into the concept of combinatorial landscapes Reidys and Stadler (2002); Schiavinotto and Stützle (2007); Tayarani-N and Prugel-Bennett (2014), and problem solving can be cast to a search over a space of states. Such problems often are very hard Cheeseman et al. (1991). In order to pursue an efficient search, it is vital to develop techniques to identify problem features and of exploiting the local search structure Frank et al. (1997); Hoffmann (2001); Hoos and Stützle (2004). In particular, it is important to reduce and decompose the problem preserving structural features that enable heuristic search and with the effective problem size factored out Slaney and Walsh (2001); Mears and de la Banda (2015). It is also significant to tackle ruggedness Billinger et al. (2014) and neutrality (plateaus) Collins (2006); Benton et al. (2010) in the landscapes.

In number theory, a $k$-tuple $\mathcal{H}_{k}$ is admissible if $\phi_{p}(\mathcal{H}_{k})<p$ for every prime $p$, where $\mathcal{H}_{k}=(h_{1},\ldots,h_{k})$ is a strictly increasing sequence of integers, and $\phi_{p}(\mathcal{H})$ denotes the number of distinct residue classes modulo $p$ occupied by the elements in $\mathcal{H}_{k}$ Goldston et al. (2009). The objective is to minimize the diameter of $\mathcal{H}_{k}$, i.e., $d(\mathcal{H}_{k})=h_{k}-h_{1}$, for a given $k$.

The early work Hensley and Richards (1974); Gordon and Rodemich (1998); Clark and Jarvis (2001) to compute narrow admissible tuples has been motivated by the incompatibility of the two long-standing Hardy-Littlewood conjectures.

Admissible sets have been used in the recent breakthrough work to find small gaps between primes. In Goldston et al. (2009), it was proved that any admissible $\mathcal{H}_{k_{0}}$ contains at least two primes infinitely often, if $k_{0}$ satisfies some arithmetic properties. In Zhang (2014), it was proved that a finite bound holds at $k_{0}\geq 3.5\times 10^{6}$. The bound was then quickly reduced to $d^{*}(\mathcal{H}_{105})=600$ Maynard (2015) and $d^{*}(\mathcal{H}_{50})=246$ Polymath (2014b), and wider ranges of $k_{0}$ also were obtained on bounded intervals containing many primes Polymath (2014b). Moreover, admissible sets have been used to find large gaps between primes Ford et al. (2015).

Most of the existing techniques to find narrow admissible tuples are sieve methods Hensley and Richards (1974); Clark and Jarvis (2001); Gordon and Rodemich (1998); Polymath (2014a, b), although a few local optimizations were proposed recently Polymath (2014a, b).

In this paper, we formally model this problem into a combinatorial optimization problem, and design search strategies to tackle the landscape, by utilizing the local search structure. Our solver is systematically tested to show its effectiveness.

2 Search Problem Formulation

For a given $k$, a candidate number set $\mathcal{V}$ with $|\mathcal{V}|\geq k$ can be precomputed, and a required prime set $\mathcal{P}$, where each prime $p\leq k$, can be determined. Each $\mathcal{H}$ is obtained by selecting the numbers from $\mathcal{V}$, and the admissibility is tested using $\mathcal{P}$.

Definition 1** (Constraint Optimization Model).**

For a given $k$, and given the required $\mathcal{V}$ and $\mathcal{P}$, the objective is to find a number set $\mathcal{H}\subseteq\mathcal{V}$ with the minimal $d(\mathcal{H})$ value, subject to the constraints $|\mathcal{H}|=k$ and $\mathcal{H}$ is admissible.

For convenience, $\mathcal{V}$, $\mathcal{P}$, and $\mathcal{H}$ are assumed to be sorted in increasing order. We denote $\mathcal{H}$ as $\mathcal{H}_{k}$ if $|\mathcal{H}|=k$, as $\tilde{\mathcal{H}}$ if it is admissible, and as $\tilde{\mathcal{H}}_{k}$ if it satisfies both of the constraints.

Given $\mathcal{V}$ and $\mathcal{P}=(p_{1},\cdots,p_{i},\cdots,p_{|\mathcal{P}|})$, the following three data structures are defined for facilitating the search:

Definition 2** (Residue Array $\mathcal{R}_{v}$).**

For $v\in\mathcal{V}$, $\mathcal{R}_{v}$ is calculated on $\mathcal{P}$: its $i$th row is $r_{v,i}=v\operatorname{mod}p_{i}$ for $i\in[1,|\mathcal{P}|]$.

Definition 3** (Occupancy Matrix $\mathcal{M}$).**

$\mathcal{M}$* is an irregular matrix, in which each row $i$ contains $p_{i}$ elements corresponding to the residue classes modulo $p_{i}$. For any given number set $\mathcal{H}\subseteq\mathcal{V}$, there is $m_{i,j}=\sum_{v\in\mathcal{H}}\mathbbm{1}(j\equiv r_{v,i}+1)$ for $j\in[1,p_{i}]$, which means the count of numbers in $\mathcal{H}$ occupying each residue class modulo $p_{i}$.*

Definition 4** (Count Array $\mathcal{F}$).**

$\mathcal{F}$* is an array, in which each row $i$ gives the count of zero elements in the $i$th row of $\mathcal{M}$.*

Property 1**.**

The space requirements for $\{\mathcal{R}_{v}|v\in\mathcal{V}\}$, $\mathcal{M}$, $\mathcal{F}$ are respectively $|\mathcal{V}|\cdot|\mathcal{P}|$, $\sum_{i\in[1,|\mathcal{P}|]}p_{i}$, and $|\mathcal{P}|$.

Figure 1 gives the example for an admissible set $\tilde{\mathcal{H}}_{7}=(0,2,8,12,14,18,30)$ with $d(\tilde{\mathcal{H}}_{7})=30$, the full prime set $\mathcal{P}$, $\mathcal{R}_{v}$ for each $v\in\tilde{\mathcal{H}}_{7}$, and the corresponding $\mathcal{M}$ and $\mathcal{F}$.

$\mathcal{H}$ and its corresponding $\mathcal{M}$ and $\mathcal{F}$ have a few properties:

Property 2** (Admissibility).**

$\mathcal{H}$* is admissible if $f_{i}>0$, $\forall i$.*

There is a constraint violation at row $i$ if $f_{i}=p_{i}-\phi_{p_{i}}\equiv 0$. The total violation count should be 0 for each $\tilde{\mathcal{H}}$.

Property 3**.**

Let $W_{i,j}=\{v\in\mathcal{H}|r_{v,i}\equiv j-1\}$, it contains all numbers occupying $(i,j)$ of $\mathcal{M}$, and there is $|W_{i,j}|=m_{i,j}$.

Property 4**.**

For each row $i$, there is $\sum_{j\in[1,p_{i}]}m_{i,j}=|\mathcal{H}|$.

There are two basic properties based on an admissible $\tilde{\mathcal{H}}_{k}$:

Property 5** (Offsetting).**

For any $c\in\mathbb{Z}$, $\mathcal{H}_{k}^{[c]}=(h_{1}+c,\ldots,h_{k}+c)$ is admissible, and there is $d(\mathcal{H}_{k}^{[c]})=d(\tilde{\mathcal{H}}_{k})$.

Property 6** (Subsetting).**

Any subset of $\tilde{\mathcal{H}}$ is admissible.

Properties 5 and 6 were observed in Polymath (2014b). Offsetting can be seen as rotating of residue classes at each row in $\mathcal{M}$; and subsetting does not decrease each row of $\mathcal{F}$.

Defining a compact $\mathcal{V}$ is nontrivial for reducing the size of problem space, which is exponential in $|\mathcal{V}|$.

One plausible way is to let $\mathcal{V}$ include all numbers in $[0,U]$, and set $h_{1}=0$. Let $d_{k}^{LB}$ and $d_{k}^{UB}$ be the best-so-far lower and upper bounds of the optimal value of $d(\tilde{\mathcal{H}}_{k})$. During the search, $|\mathcal{V}|=U$ can be bounded by $d_{k}^{UB}$. However, it appears that $d_{k}^{UB}$ is very close to $\lfloor k\log k+k\rfloor$ Polymath (2014b), which might still be very large when $k$ is big.

Based on Property 4, as $p_{i}$ is small, the average $m_{i,j}$ would be large. For the rows with $f_{i}=1$, it would be difficult to find a useful heuristic for changing the unoccupied column $j$ at row $i$. Intuitively, each unoccupied location $(i,j)$ in $\mathcal{M}$ can be assumed to be unchanged during the search. Thus, sieving can be applied to remove any numbers in $[0,U]$ that occupy those unoccupied locations, which could be found using Property 3, for a set of the smallest primes $\mathcal{P}_{R}\in\mathcal{P}$. After sieving, the proportion of remaining numbers is $\alpha\approx 1-\prod_{p_{i}\in\mathcal{P}_{R}}(1-1/p_{i})$. The completeness of the problem space on the other combinations in $\mathcal{M}$ can be retained using a sufficiently large $U$, based on the principle behind Property 5, i.e., offsetting as a choice of residue classes.

Remark 1**.**

The original problem can be converted into a list of subproblems, where each subproblem takes each $v\in[0,U-d_{k}^{LB}]\cap\mathcal{V}$ as the starting point $h_{1}$ to obtain the minimal diameter for each $\tilde{\mathcal{H}}_{k}\subseteq\mathcal{V}_{v}=[v,v+d_{k}^{UB}]\cap\mathcal{V}$. The optimal solution is then the best solution among all subproblems.

Decomposition Friesen and Domingos (2015) has been successfully used in AI for solving discrete problems. The new perspective of the problem has two features. First, each subproblem has a much smaller state space, as $|\mathcal{V}_{v}|\approx\alpha\cdot d_{k}^{UB}$. Second, for totally around $\alpha\cdot(U-d_{k}^{LB})$ subproblems, the good solutions of neighboring subproblems might share a large proportion of elements, providing a very useful heuristic clue for efficient adaptive search in this dimension.

Let $\mathcal{P_{C}}$ contain all primes $p\leq k$. In theory, $\mathcal{P}=\mathcal{P_{C}}$, but we can reduce it to an effective subset. Based on Property 4, if $p_{i}$ is large, the average $m_{i,j}$ would be small. Some rows of $\mathcal{F}$ would always have $f_{i}>0$, even for $\mathcal{H}=\mathcal{V}$. The set of corresponding primes, named $\mathcal{P}_{L}$, thus can be removed from $\mathcal{P}$, without any loss of the completeness for testing the admissibility. The effective prime set would be $\mathcal{P}=\mathcal{P_{C}}-\mathcal{P}_{L}$.

Algorithm 1 gives the specific realization for obtaining $\mathcal{V}$ and $\mathcal{P}$. Here $\mathcal{P}_{R}$ in Line 1 and $\mathcal{V}$ in Line 2 are obtained using the setting in the greedy sieving method Polymath (2014b).

3 Search Algorithm

In this section, the basic operations on auxiliary data structures are first introduced. Some search operators are then realized. Finally, we describe the overall search algorithm.

3.1 Operations on $\mathcal{R}_{v}$, $\mathcal{M}$ and $\mathcal{F}$

For every $v\in\mathcal{V}$, $\mathcal{R}_{v}$ is calculated in advance. For $\mathcal{H}\subseteq\mathcal{V}$, the corresponding $\mathcal{M}$ and $\mathcal{F}$ are synchronously updated locally. The space requirements are given in Property 1.

There are two elemental 1-move operators, i.e., adding $v\notin\mathcal{H}$ into $\mathcal{H}$ to obtain $\mathcal{H}=\mathcal{H}+\{v\}$, and removing $v\in\mathcal{H}$ from $\mathcal{H}$ to obtain $\mathcal{H}=\mathcal{H}-\{v\}$, for given $\mathcal{H}\subseteq\mathcal{V}$ and $v\in\mathcal{V}$.

Property 7** (Connectivity).**

The two elemental 1-move operators possess the connectivity property for each $\mathcal{H}\in\mathcal{V}$.

The connectivity property Nowicki and Smutnicki (1996) states that there exists a finite sequence of such moves to achieve the optimum state from any state in the search space.

For $\mathcal{H}\equiv\varnothing$, there are $m_{i,j}=0$ for each $i,j$, $f_{i}=p_{i}$ for each $i$, based on Definitions 3 and 4. The $\mathcal{M}$ and $\mathcal{F}$ for any $\mathcal{H}$ can be constructed by adding each $v\in\mathcal{H}$ using Algorithm 2. For any two states $\mathcal{H}_{A}$ and $\mathcal{H}_{B}$, $\mathcal{H}_{A}$ can be changed into $\mathcal{H}_{B}$ by adding each $v\in\mathcal{H}_{B}-\mathcal{H}_{A}$ and by removing each $v\in\mathcal{H}_{A}-\mathcal{H}_{B}$. The total number of 1-moves is $L=|\mathcal{H}_{A}\cup\mathcal{H}_{B}|-|\mathcal{H}_{A}\cap\mathcal{H}_{B}|$, i.e., which can be seen as the distance Reidys and Stadler (2002) between two states. The shorter the distance, the more similar the two states are.

Algorithms 2 and 3 respectively give the operations of updating $\mathcal{M}$ and $\mathcal{F}$ for the two elemental 1-move operators.

Property 8** (Time Complexity).**

Algorithms 2 and 3 have the time complexity $O(|\mathcal{P}|)$ in updating $\mathcal{M}$ and $\mathcal{F}$.

In the following realizations, we will focus on the search transitions between $\tilde{\mathcal{H}}$ states. The admissibility testing (Property 2) on each $\tilde{\mathcal{H}}$ is not explicitly applied. Instead, VioCheck in Algorithm 4 is used to check $\Delta$, i.e., the violation count to be increased, if adding $v$ into $\mathcal{H}$, using $R_{v}$, $\mathcal{M}$ and $\mathcal{F}$.

3.2 Search Operators

We first realize some elemental and advanced search operators to provide the transitions between admissible states.

3.2.1 Side Operators

Let $Side$={Left, Right} define the two mutually reverse sides of ${\mathcal{H}}$. For an admissible state $\tilde{\mathcal{H}}$, each side operator tries to add or remove a number at the given side of $\tilde{\mathcal{H}}$ to obtain the admissible tuple with a diameter as narrow as possible.

SideRemove just removes the element at the given $Side$ from $\tilde{\mathcal{H}}$, and its output is admissible, according to Property 6.

Algorithm 5 defines the operation SideAdd for adding a number $v$ into $\tilde{\mathcal{H}}$. To retain the admissibility, the number to be added is tested using Algorithm 4 to ensure the admissibility.

3.2.2 Repair Operator

The Repair operator repairs $\tilde{\mathcal{H}}$ into $\tilde{\mathcal{H}}_{k}$ using side operators: While $|\tilde{\mathcal{H}}|<k$, the SideAdd operator is iteratively applied on each side of $\tilde{\mathcal{H}}$, and the better one is kept; While $|\tilde{\mathcal{H}}|>k$, the SideRemove operator is iteratively applied on each side of $\tilde{\mathcal{H}}$, and the better one is kept. Finally, $\tilde{\mathcal{H}}_{k}$ is obtained as $|\tilde{\mathcal{H}}|=k$.

3.2.3 Shift Search

Algorithm 6 gives the realization of the ShiftSearch operator. The side is selected at random (Line 2). Each shift Polymath (2014a) is realized by combining SideRemove and SideAdd (Line 4), leading to a distance of 2 to the original state. Starting from $\tilde{\mathcal{H}}_{O}$, we applies totally up to $N_{L}$ shifts (Line 3) unless SideAdd fails (Line 5), and the best state is kept as $\tilde{\mathcal{H}}_{N}$ (Line 6). The state $\tilde{\mathcal{H}}_{N}$ is accepted immediately if $d_{N}\leq d_{O}$, or with a probability otherwise (Line 8), following the same principle as in simulated annealing Kirkpatrick et al. (1983).

3.2.4 Insert Moves

Algorithm 7 gives the realization of the InsertMove operator to work on the input $\tilde{\mathcal{H}}$ for obtaining an admissible output.

The operator is realized in three levels, as defined by the parameter $Level\in\{0,1,2\}$. For each $v$ in a compact set $\mathcal{V}_{in}=[h_{1},h_{|\tilde{\mathcal{H}}|}]\cap\mathcal{V}-\tilde{\mathcal{H}}$ (Line 2), the violation count $\Delta$ is calculated using Algorithm 4 (Line 3). At level 0, the value $v$ is immediately inserted into $\tilde{\mathcal{H}}$ if $\Delta\equiv 0$ (Line 4). Otherwise, if $Level>0$ and $\Delta\equiv 1$, the violation row $i$ is found using VioRow (Line 5), which is simply realized by returning $i$ as the conditions are satisfied at Line 3 of Algorithm 4, and then $v$ is stored into the set $Q_{i}$ (Line 5) starting from $\varnothing$ (Line 1).

At levels 1 and 2, we compare $|Q_{i}|$ and $m_{i,sb}$, where $sb$ is the second best location in row $i$ of $\mathcal{M}$. Based on Property 3, $|W_{i,sb}|=m_{i,sb}$. Note that the admissibility is retained after adding elements in $Q_{i}$ and removing elements in $W_{i,sb}$.

Remark 2**.**

For $\tilde{\mathcal{H}}$=InsertMove($\tilde{\mathcal{H}}$), $d(\tilde{\mathcal{H}})$ is non-increasing at all levels. $|\tilde{\mathcal{H}}|$ is respectively increased by 1 and $|Q_{i}|-m_{i,sb}$ at levels 0 and 1, and keeps unchanged at level 2.

In general, InsertMove is successful if it can increase $|\tilde{\mathcal{H}}|$. However, the neighborhood might contains too many infeasible moves, as many 1-moves would trigger multiple violations. It might be inefficient to use systematic adjustments Polymath (2014a). Our implementation targets on feasible moves intelligently by utilizing the violation check clues.

3.2.5 Local Search

Algorithm 8 gives the realization of the LocalSearch operator. We will only focus on the case of improving the input state with $|\tilde{\mathcal{H}}|=k$. Let the input have $d_{0}=d(\tilde{\mathcal{H}})$. The SideRemove operator is first applied for $N_{S}$ times (Lines 1-3). Its output has $|\tilde{\mathcal{H}}|<k$, and $d_{1}=d(\tilde{\mathcal{H}})<d_{0}$. The InsertMove operator is then applied for up to $N_{I}$ times (Lines 4-6). Based on Remark 1, the output has $d_{2}=d(\tilde{\mathcal{H}})\leq d_{1}$. If this step leads to $|\tilde{\mathcal{H}}|\geq k$, the final output after repairing definitely has a lower diameter than $d_{0}$. Otherwise, it is still possible to produce a better output as the state is being repaired (Line 7).

3.3 Region-based Adaptive Local Search (RALS)

The region-based adaptive local search (RALS) is realized to tackle the problem decomposition as described in Remark 1.

Let us consider the problem along the dimension of the numbers in $\mathcal{V}$. For each $\tilde{\mathcal{H}}_{k}$, it can be indexed by $(h_{1},d(\tilde{\mathcal{H}}_{k}))$. Let $f^{*}({v})$ be the optimal diameter for all $\tilde{\mathcal{H}}_{k}$ with $h_{1}=v$, we can form a set of points $\{(v_{1},f^{*}(v_{1})),\cdots,(v_{|\mathcal{V}|},f^{*}(v_{|\mathcal{V}|})\}$. It can be seen as a one-dimensional fitness landscape representing the fitness function $f^{*}(v)$ on the discrete variable from $v\in\mathcal{V}$. Note that the optimal solution on this fitness landscape is the optimal solution of the original problem.

Essentially, we would like to focus the search on those promising regions where $f^{*}(v)$ has higher quality. Nevertheless, early search can provide some clues for narrowing down promising regions, even though the fitness landscape itself is not explicit at the beginning, as $f^{*}(v)$ at each $v$ can be revealed through extensive local search.

3.3.1 Database Management

We use a simple database, denoted by DB, to keep the high-quality solutions $\tilde{\mathcal{H}}_{k}$ found during the search, and index each of them as $(v,f({v}))$, where $v=h_{1}$, and $f({v})=d(\tilde{\mathcal{H}}_{k})$ for each $\tilde{\mathcal{H}}_{k}$. For each $v$, only the best-so-far solution and the corresponding $f({v})$ is kept. Here $f({v})$ plays the role of a virtual fitness function that is updated during the search process to approximate the real fitness function $f^{*}(v)$.

The database is managed in a region-based mode. Specifically, the total range $[0,U]$ of the numbers in $\mathcal{V}$ is divided into $N_{R}$ regions. There are three basic operations for DB.

The DBInit operator is used for providing the initialization. The greedy sieve Polymath (2014b) is applied to generate a state $\tilde{\mathcal{H}}_{k}$ in each region for forming the initial $f({v})$.

The DBSelect operator is used for selecting one incumbent state to apply the search operation. In the region-based mode, there are two steps to provide the selection. In the first step, each region provides one candidate. In this paper, we greedily choose the best solution in each region. In the second step, the incumbent state is selected from the candidates provided by all regions. We consider the following implementation. At the probability $\gamma$, the candidate is selected at random. Otherwise, tournament selection is applied to select the best solution among totally $N_{T}$ randomly chosen candidates.

The DBSave operator simply stores each $\tilde{\mathcal{H}}_{k}$ into DB, and updates $f({v})$ internally. Dominated solutions are discarded.

3.3.2 Algorithm Realization

Algorithm 9 gives the implementation of RALS to obtain $\mathcal{H}_{k}^{*}$ for a given $k$. First, $\mathcal{V}$ and $\mathcal{P}$ are initialized using Algorithm 1, using $U=\lceil 1.5\cdot(k\log k+k)\rceil$ for the range $[0,U]$ (Line 1) to ensure $U>d_{k}^{UB}$. Afterward, the database DB with $N_{R}$ regions is initiated using the DBInit operator (Line 2).

The search process runs $T$ iterations in total. In each iteration, we first select one incumbent solution $\tilde{\mathcal{H}}_{k}$ from DB using the DBSelect operator (Line 4). Then the actual search tackles two parts of the problem. The ShiftSearch operator is used to search on the virtual fitness landscape $f({v})$ (Line 5). The LocalSearch operator is then applied to improve $f({v})$ locally (Lines 6-7). For each search operator in Lines 5-7, the DBSave operator is applied to store newly generated solutions. Finally, the best solution $\mathcal{H}^{*}$ in DB is returned.

4 Results and Discussion

We now turn to the empirical evaluation of the proposed algorithm. For the benchmark instances, we refer to an online database Sutherland (2015) that has been established and extensively updated to contain the narrowest admissible $k$-tuples known for all $k\leq 5000$. The algorithm is coded in Java, and our experiments were run on AMD 4.0GHz CPU. For each instance, 100 independent runs were performed.

4.1 Results by Existing Methods

Most existing techniques to solve this problem are constructive and sieve methods Polymath (2014a). The sieve methods are realized by sieving an integer interval of residue classes modulo primes $p<k$ and then selecting an admissible $k$-tuple from the survivors. The easiest way to construct a narrow $\tilde{\mathcal{H}}_{k}$ is using the first $k$ primes past $k$ Zhang (2014). As an optimization, the sieve of Eratosthenes takes $k$ consecutive primes, to search starting from $p<k$, in order to select one among the admissible tuples that minimize the diameter.

The Hensley-Richards sieve Hensley and Richards (1974) uses a heuristic algorithm to sieve the interval $[-x/2,x/2]$ to obtain $\tilde{\mathcal{H}}_{k}$, leading to the upper bound Polymath (2014b):

[TABLE]

The Schinzel sieve, as also considered in Gordon and Rodemich (1998); Clark and Jarvis (2001), sieves the odd rather than even numbers. In the shifted version Polymath (2014b), it sieves an interval $[s,s+x]$ of odd integers and multiples of odd primes $p\leq p_{m}$, where $x$ is sufficient large to ensure at least $k$ survivors, and $m$ is sufficient large to ensure that the survivors form $\tilde{\mathcal{H}}$, $s\in[-x/2,x/2]$ is the starting point to choose for yielding the smallest final diameter.

As a further optimization, the shifted greedy sieve Polymath (2014b) begins as in the shifted Schinzel sieve, but the minimally occupied residue class are greedily chosen to sieve for primes $p>\tau\sqrt{k\log k}$, where $\tau$ is a constant. Empirically, it appears to achieve the bound Polymath (2014a):

[TABLE]

Table 1 lists the upper bounds obtained by applying a set of existing techniques, including $k$ primes past $k$, Eratosthenes (Zhang) sieve, Hensley-Richards sieve, Schinzel and Shifted Schinzel sieve, by running the code111http://math.mit.edu/~drew/ompadm_v0.5.tar provided in Polymath (2014b), on $k=\{1000,2000,3000,4000,5000\}$. The best known results are retrieved from Sutherland (2015).

4.2 Results by RALS algorithm

Table  2 lists the average results of different versions of the proposed RALS algorithm. The “BaseVer” version is defined with the following settings. For the database DB, we use $N_{R}=20$. For its DBSelect operator, there are $\gamma=0.01$ and $N_{T}=4$. For the search loop, we consider $T=1000$ iterations. For the ShiftSearch operator, we set $N_{L}=10$ and ${\beta}=1$. For the InsertMove operator, there is $Level=2$. For the parameters of LocalSearch in Algorithm 9, we set $N_{I1}=500$ and $N_{I2}=0$. The other versions are then simply the “BaseVer” version with different parameters.

With $T=0$, the algorithm returns the best results obtained by the shifted greedy sieve Polymath (2014a) in the $N_{R}$ regions. The results are significantly better than the sieve methods in Table 1. The search operators in RALS show their effectiveness as all RALS versions with $T>0$ perform significantly better than the version with $T=0$.

Note that “BaseVer” has $Level=2$, we can compare the RALS versions with different levels $\{0,1,2\}$ in the InsertMove operator of LocalSearch. On the performance, the version with a higher level produces better results than that of a lower level. With greedy search only, the first LocalSearch works as an efficient contraction process Polymath (2014a). As described in Remark 2, InsertMove performs greedy search at levels 0 and 1, but performs plateau moves at level 2, from the perspective of updating $|\tilde{\mathcal{H}}|$. At level 0, the search performs elemental 1-moves. At level 1, the search can be in a very large neighborhood although it has a low time complexity. Plateau moves is used at level 2 to find exits, as remaining feasible moves are more difficult to check. Finding exits to leave plateaus Hoffmann (2001); Frank et al. (1997) has been an important research topic about the local search topology on many combinatorial problems Bonet and Geffner (2001); Benton et al. (2010); Sutton et al. (2010). From the viewpoint of the LocalSearch operator, the plateau moves on the part solved by InsertMove help escaping from local mimina in the landscape of the subproblem.

We compare the RALS versions with different $N_{I1}\in\{100,500,1000\}$, as “BaseVer” has $N_{I1}=500$. Especially for $k\in\{3000,4000\}$, the improvements of using higher $N_{I1}$ are extremely significant as $N_{I1}$ increases from 100 to 500, but are less significant as $N_{I1}$ further increases to 1000.

In “BaseVer”, the second LocalSearch in Line 7 of Algorithm 9 is actually not used if it has $N_{I2}=0$. As we increase $N_{I2}$ to $10$, the instance $k=1000$ can be fully solved to the best known solution, and the instance $k=5000$ can also be solved to obtain a significantly better result.

Lines 1-3 of Algorithm 8 might be viewed as perturbation, an effective operator in stochastic local search Hoos and Stützle (2004) to escape from local minima on rugged landscapes Tayarani-N and Prugel-Bennett (2014); Billinger et al. (2014). In RALS, the second LocalSearch essentially applies a larger perturbation than the first LocalSearch.

Table 2 also gives the comparison for RALS with $\gamma\in\{0.001,0.01,0.1,1\}$ for selecting the incumbent state in DBSelect. The larger the $\gamma$, the more random the selection is. The best performance is achieved as $\gamma=0.1$, neither too greedy nor too random. We can gain some insights from a typical snapshot of the virtual fitness landscape $f(v)$, as shown in Figure 2. It is easy to spot the valley with high-quality solutions, as they provide significant clues for adaptive search. Meanwhile, the noises on the fitness landscape might reduce the effectiveness of pure greedy search. Thus, there is a trade-off between greedy and random search.

Tables 3 lists the performance measures, including the average, the successful rate of finding best known solutions (SuccRate), and the calculating time, by the “BaseVer” version with both $\gamma=0.1$ and $N_{I2}=10$, as $T=100$ and $T=1000$. This version achieves high SuccRate for $k\in\{1000,2000,3000\}$, and moderate SuccRate for $k\in\{4000,5000\}$, as $T=1000$. It also reaches reasonable good SuccRate as $T=100$, with a lower execution time.

Finally, we apply RLAS to compare the results for $k\in[2500,5000]$ in Sutherland (2015). In Table 4, we list the new upper bound $\mathcal{H}_{k}^{*}$ and the improvement on the diameter $\delta d$ for each $k$ of the 48 instances. Eight instances among them have $\delta d\geq 10$. Thus, AI-based methods might make further contributions to pure mathematics applications.

5 Conclusions

We presented a region-based adaptive local search (RALS) method to solve a case of pure mathematics applications for finding narrow admissible tuples. We formulated the original problem into a combinatorial optimization problem. We showed how to exploit the local search structure to tackle the combinatorial landscape, and then to realize search strategies for adaptive search and for effective approaching to high-quality solutions. Experimental results demonstrated that the method can efficiently find best known and new solutions.

There are several aspects of this work that warrant further study. A deeper analysis might be applied to better identify properties of the local search topology on the landscape. One might also apply advanced AI strategies, e.g., algorithm portfolios Gomes and Selman (2001) and SMAC Hutter et al. (2011), to obtain an even greater computational advantage.