Dynamic thresholding search for the feedback vertex set problem

Basic definitions

Given a directed graph $G=(V,E)$, basic definitions that are useful for describing the proposed IDTS algorithm are presented as below.

Definition 1: a $criticalvertex$ of G is a vertex that belongs to a FVS. We use C to denote the set of critical vertices that have been detected. C is a FVS only when all vertices of the FVS are detected.

Definition 2: an $uncriticalvertex$ is a vertex that does not belong to a FVS. We use U to denote the set of uncritical vertices, and $V=C\cup U,C\cap U=\varnothing$.

Definition 3: a $redundantvertex$ refers to a vertex that is recognized as critical or uncritical according to the exact rules proposed by Levy & Low (1988). We use $V_r$ to denote the set of redundant vertices that have been detected, $C_r$ to denote the set of critical vertices of $V_r$, $U_r$ to denote the set of uncritical vertices of $V_r$, and $V_r=C_r\cup U_r,C_r\cap U_r=\varnothing$.

Definition 4: $V_0$ refers to the set of residual vertices after applying the removal exact algorithm proposed by Levy & Low (1988). $C_0$ denotes the set of feedback vertices of $V_0$ (that is, all vertices of a FVS are detected and belong to $C_0$), $U_0$ denotes the set of non-feedback vertices of $V_0$, and $V_0=C_0\cup U_0,C_0\cap U_0=\varnothing$. Levy & Low (1988) proved that the FVS of the reduced graph plus the FVS removed in the reduction process composes the FVS of the original graph. Let $C_r$ be the set of vertices removed in the reduction process, and $C_0^{\ast }$ be the minimum FVS of the reduced graph $G=(V_0,E_0)$, where $E_0=V_0\times V_0\cap E$. Then, $C_r\cap C_0^{\ast }=\varnothing$ and $C_0^{\ast }\cup C_r$ is a minimum FVS of the given graph $G=(V,E)$. In this case, only the feedback vertices of the reduced graph need to be found out.

In summary, the vertex set V of G consists of two disjoint sets $\{V_0,V_r\}$ or four disjoint sets $\{C_0,U_0,C_r,U_r\}$.

Definition 5: a directed acyclic graph (DAG) (Bangjensen & Gutin, 2008) is a directed graph with no directed cycles.

Vertices in each DAG are in a topological ordering where the starting point of every directed edge is ahead of its terminal point (Galinier, Lemamou & Bouzidi, 2013). For a vertex set $U\subset V$, we notice that the induced subgraph $G_U=(U,E_U)$, $E_U=U\times U\cap E$ is acyclic if and only if $V\mathrm{\setminus }U$ is a FVS. Hence, the objective of the FVSP is to find the set U that has the maximum cardinality to make $G_U=(U,E_U)$ acyclic.

Figure 2 presents an example illustrating these basic definitions. For the given graph $G=(V,E)$, let $\{a,b,d,h\}$ be the current FVS. The set of critical vertices C is the current FVS $\{a,b,d,h\}$, and the set of remaining vertices is the set of uncritical vertices U, i.e., $\{c,e,f,g,i,j\}$. According to the rules in “Reduction Procedure”, $h$ can be recognized as a critical vertex ( $C_r=\{h\}$, purple vertex) and $f,g,i,j$ as uncritical vertices ( $U_r=\{f,g,i,j\}$, dark blue vertices). Thus the set of redundant vertices $V_r$ is $\{f,g,h,i,j\}$, and the set of residual vertices $V_0$ is $\{a,b,c,d,e\}$. $V_0$ can be divided into $C_0=\{a,b,d\}$ (orange vertices) and $U_0=\{c,e\}$ (blue vertices). Clearly, the graph induced by the vertices in $U=\{c,e,f,g,i,j\}$ is a DAG without directed cycles.

Solution representation and fitness function

The solution representation and fitness function of the FVSP are given as follows.

Solution representation: the constraint of the FVSP is that there is no cycle in $U_0$ after removing the set of redundant vertices $C_r\cup U_r$ and the set of critical vertices $C_0$. To quickly assess the number of cycles in $U_0$ after each neighborhood operation, the number of conflicts (see “Preliminaries”) is taken as the number of cycles (Galinier, Lemamou & Bouzidi, 2013). Let $\pi$ be an assignment of the vertices of $U_0$ to the positions $\{1,2,\dots ,|U_0|\}$, and thus the permutation $\pi$ denotes the candidate solution (Galinier, Lemamou & Bouzidi, 2013).

Fitness function: To evaluate the quality of the FVS C, the evaluation or fitness function counts the number of vertices in C. Recall that $U_r$ is the set of uncritical vertices of the redundant vertices and $\pi$ is the corresponding permutation solution of C. The fitness function $f_0$ (to be minimized) is given by

(1) $$Minimum{f}_0(\pi )=|V|-|U_r|-|\pi |$$

Thus, the minimization of the function $f_0$ is equal to the maximization of the fitness function $f$, which is expressed as:

(2) $$Maximizef(\pi )=|\pi |$$

Basic steps

This section introduces the iterated dynamic thresholding algorithm for solving the FVSP, which is composed of five main procedures as shown in Algorithm 1.

Reduction procedure: IDTS adopts a set of conventional reduction rules (Levy & Low, 1988) to simplify the given graph G. Firstly, a set of redundant vertices $V_r$ (made up of the set of critical vertices $C_r$ and the set of uncritical vertices $U_r$) is confirmed according to those rules. Then, the redundant vertices and related edges (whose starting or ending vertex is a redundant vertex) are deleted, reducing the input graph $G=(V,E)$ to the reduced graph $G=(V_0,E_0)$ (see “Reduction Procedure”).

Initialization procedure: this procedure greedily chooses a vertex (Cai, Huang & Jian, 2006) such that its insertion to $\pi$ does not increase the number of cycles. This process continues until no such vertex can be inserted (see “Greedy Initialization”).

Local search procedure: this procedure consists of two complementary search stages: a dynamic thresholding search stage (diversification) to extend the search to unexplored regions, and a descent search stage (intensification) to find new local optimal solutions with improved quality. These two stages alternate until the best-found solution cannot be further improved for $\omega$ continuous local search rounds (see “Local Search”).

Perturbation procedure: When the search is considered as trapped in a deep local optimum, the perturbation procedure is initiated to move some specifically identified vertices between $U_0$ and $C_0$ to relieve the search from the trap. The solution perturbed is then adopted to start the next round of the local search procedure (see “Perturbation Procedure”).

Recovery procedure: If the best solution ever found cannot be improved after $\gamma$ continuous local search rounds and the perturbation phase, the search then terminates and the recovery procedure starts. The best solution (minimum feedback vertex set) found in the search procedure is recorded as $C_0^{\ast }$ and returned as the input of the recovery procedure. The current reduced graph $G=(V_0,E_0)$ is recovered to the original graph $G=(V,E)$, and the FVS $C_0^{\ast }$ for $G=(V_0,E_0)$ is correspondingly projected back to $C_0^{\ast }\cup C_r$ for $G=(V,E)$ (see “Recovery Procedure”).

Reduction procedure

The graph reduction procedure follows three rules when traversing all vertices in a given graph $G=(V,E)$ and processes those that satisfy any rule proposed by Levy & Low (1988). The three reduction rules are as follows.

Rule 1: If the in-degree (out-degree) of a vertex $v$ is 0, that is, $v$ is an uncritical vertex, then $v$ and all its edges can be deleted without missing any optimal feedback vertex of G. Such vertices are added into the redundant uncritical set $U_r$ ( $U_r=U_r\cup \{v\}$). For example, as shown in Fig. 4B, the edges of the vertex $g$ whose in-degree is 0, and those of the vertex $j$ whose out-degree is 0 can be deleted.

Rule 2: If the in-degree (out-degree) of a vertex $v$ is 1, and there is no self-loop, then vertex $v$ can be merged with the unique precursor (successor) vertex without missing any optimal feedback vertex of G. The merging process is that all edges connected to the vertex $v$ are linked to the unique precursor (successor) vertex of $v$. Such vertices are added into the redundant uncritical set $U_r$ ( $U_r=U_r\cup \{v\}$). Figure 4C displays that the vertex $f$ with an in-degree of 1 can be merged with the precursor vertex $a$, and the vertex $i$ with an out-degree of 1 can be merged with the successor vertex $b$.

Rule 3: If a self-loop exists for a vertex $v$, then $v$ and all its edges can be deleted and recovered as a part of the feedback vertex set without losing any optimal feedback vertex of G. Such vertices are added into the redundant uncritical set $C_r$ ( $C_r=C_r\cup \{v\}$). As shown in Fig. 4D, the self-loop vertex $h$ and its connected edges can be deleted.

After deleting the sets $U_r$ and $C_r$, the remaining vertex set $V_0=V\mathrm{\setminus }(U_r\cup C_r)$, and the reduced subgraph $G=(V_0,E_0)$, $(E_0=(V_0\times V_0\cap E))$. After obtaining reduced G, the greedy initialization is used to generate an initial solution for it.

Greedy initialization

Given the reduced subgraph $G=(V_0,E_0)$, its critical vertex set is defined as $C_0$, and the uncritical vertex set as $U_0$. Recall that $\pi$ is an assignment of the vertices of $U_0$ to the positions $\{1,2,\dots ,|U_0|\}$. We initialize $C_0=V_0$, $\pi =\mathrm{\varnothing }$. Then, $\pi$ is iteratively extended in the greedy procedure by inserting a minimum-score vertex $v$ of $C_0$ until no vertex can be inserted.

(1) Calculate $N_I(\pi )$: recall that $N_I(\pi )$ contains the vertices that satisfy $g(\pi \oplus <v,i_{-}(v),i_{+}(v)>)=0$. Thus, we only has to traverse all vertices and add the vertex whose $g(\pi \oplus <v,i_{-}(v),i_{+}(v)>)=0$ into $N_I(\pi )$.

(2) Select one vertex: choose a vertex $v$ with the minimum score (Eq. (6), Cai, Huang & Jian, 2006) in $N_I(\pi )$ and insert it into the current permutation $\pi \leftarrow \pi \oplus <v,i_{-}(v),i_{+}(v)>$.

(6) $$\begin{array}{r}score(v)=|de{g}^{-}(v)+de{g}^{+}(v)|-\lambda \times |de{g}^{-}(v)-de{g}^{+}(v)|\end{array}$$where $de{g}^{-}(v)$ and $de{g}^{+}(v)$ are the in-degree and the out-degree of $v$ respectively, and $\lambda$ is a parameter ( $\lambda =0.3$ according to Cai, Huang & Jian (2006)).

(3) Update $N_I(\pi )$: after inserting the vertex $v$ into $\pi$, we only have to recalculate the number of conflicts $g(\pi \oplus <u,i_{-}(u),i_{+}(u)>)$ of its neighbors $u\in N_I(\pi )$ according to Eq. (4). Any vertex satisfying $g(\pi \oplus <u,i_{-}(u),i_{+}(u)>)\ne 0$ will be eliminated from $N_I(\pi )$. Thus, the complexity of this updating step is $O(d_{max}^2)$, where $d_{max}$ is the largest degree of a vertex in the graph. This process continues until no vertex can be inserted into $\pi$, i.e., $N_I(\pi )=\mathrm{\varnothing }$.

In this initialization, a legal conflict-free $\pi$ with a certain quality can be obtained, and further improved in the dynamic thresholding search stage of the algorithm.

To explain this process, we consider the reduced graph $G=(V_0,E_0)$ in Fig. 4D (with vertices $\{a,b,c,d,e\}$). For Fig. 4D, $C_0=\{a,b,c,d,e\}$, $U_0=\pi =\mathrm{\varnothing }$ and $g(\pi \oplus <v,i_{-}(v),i_{+}(v)>)=0$, $\mathrm{\forall }v\in C_0$. Figure 5 shows how the greedy procedure works. Firstly we calculate $N_I(\pi )=\{a,b,c,d,e\}$. Then, it is detected that $score(a)=4.7,score(b)=4,score(c)=4,score(d)=4.7,score(e)=4$ according to Eq. (6). Finally, we select a vertex from $N_I(\pi )$ with the minimum score and insert it into $\pi$. Suppose we select vertex $e$ for insertion. $N_I(\pi )$ is updated to $\{a,c,d\}$. As shown in Fig. 5A, the solution after the first greedy insertion is the permutation $\pi =\{e\}$. By repeating the above steps until $N_I(\pi )=\mathrm{\varnothing }$, we obtain the local optimal permutation $\pi =\{e,c\}$ (Fig. 5B). Meanwhile, this process may unfortunately misclassify critical vertices in permutation $\pi$. For example, the optimal permutation of Fig. 4D is $\{d,e,a\}$ with the set of critical vertices $\{b,c\}$, while the local optimal permutation of Fig. 5B is $\{e,c\}$ with the set of critical vertices $\{a,b,d\}$.

Local search

The local optimization aims to improve the initial permutation provided by the greedy initialization, and it consists of two stages. The first stage (dynamic thresholding search) brings diversity as it accepts equivalent or worse solutions (line 4 of Algorithm 3), and the second stage applies a descent search that accepts only better solutions (line 5 of Algorithm 3) to guarantee a concentrated and directed search. These two stages alternate until the best found solution cannot be further improved for $\omega$ local search rounds.

The dynamic thresholding search stage

There are many successful applications of dynamic thresholding search (Dueck & Scheuer, 1990; Moscato & Fontanari, 1990) (e.g., the frequency assignment (Diane & Nelson, 1996), quadratic multiple knapsack problem (Chen & Hao, 2015), heterogeneous fixed fleet vehicle routing problem (Tarantilis, Kiranoudis & Vassiliadis, 2004), and others (Chen & Hao, 2019; Lai et al., 2022; Zhou, Hao & Wu, 2021)).

Algorithm 2: :

Greedy initialization for FVSP

Input: A reduced graph $G=(V_0,E_0)$

Output: A solution permutation π

1  $\pi \leftarrow \mathrm{\varnothing }$

2  $N_I(\pi )\leftarrow \mathrm{\varnothing }$

3 for each $v\in V_0$ do

4   Calculate $N_I(\pi )$ according to Eq. (5)      /*Step 1*/

5 end

6 while $N_I(\pi )\ne \mathrm{\varnothing }$ do

7   Choose a vertex $v\in N_I(\pi )$ with the minimum score according to Eq. (6) and insert it into π      /*Step 2*/

8   Update $N_I(\pi )$     /*Step 3*/

9 end

10 return π

DOI: 10.7717/peerj-cs.1245/table-11

Algorithm 3: :

Local search

Input: Reduced graph $G=(V_0,E_0)$, solution π, best solution ${\pi }^{\star }$, search depth ω

Output: Improved solution π

1  $NoImprove\leftarrow 0$      /*Indicate the times of ${\pi }^{\star }$ being improved*/

2 while $NoImprove<\omega$ do

3    $NoImprove\leftarrow NoImprove+1$

4    $(\pi ,{\pi }^{\star },NoImprove)\leftarrow Dynamic\mathrm{\_}thresholding\mathrm{\_}search(\pi ,{\pi }^{\star },NoImprove)$    /*Section 3.4.1*/

5    $(\pi ,{\pi }^{\star },NoImprove)\leftarrow Descent\mathrm{\_}search(\pi ,{\pi }^{\star },NoImprove)$          /*Section 3.4.2*/

6 end

7 return π

DOI: 10.7717/peerj-cs.1245/table-12

In this work, three basic move operators (DROP, INSERT and SWAP) are adopted in the thresholding search stage that accepts both equivalent and better solutions. DROP deletes a vertex from the current permutation $\pi$ and move it to $C_0$; INSERT extends the current permutation $\pi$ by introducing a new vertex; SWAP deletes a vertex $v$ from the current permutation $\pi$ and inserts a vertex $u$ into $\pi$.

Based on these three move operators, dynamic thresholding search (DTS) adopts both the vertex-based strategy and the prohibition mechanism to balance the exploration and exploitation of the search space. In each search round at this stage, the algorithm first randomly visits all vertices in $V_0$ one by one. For each vertex $v$ considered, the set of candidate move operators is executed depending on whether the vertex is in or out of the permutation $\pi$. If the objective value of the obtained solution is not worse than the best solution found ever to a certain threshold, the move operation on the vertex is executed. Otherwise, the operation is rejected. Each time a move is taken, the concerned vertex is marked as tabu and forbidden to be moved again during the next $tt$ iterations ( $tt$ is the tabu tenure). This process continues until all vertices of $V_0$ are traversed.

Algorithm 4: :

The dynamic thresholding search

Input: Reduced graph $G=(V_0,E_0)$, solution π, best solution ${\pi }^{\star }$, the iterations without improvement NoImprove

Output: Solution π, best solution ${\pi }^{\star }$, the iterations without improvement NoImprove

1 Randomly shuffle all vertices in V0

2 for each $v\in V_0$ do

3   if $v\notin \pi$ then

4     if $thenumberofconflictsafterinsertingvinto\pi is0$ then

5           Insert v into π      /* INSERT operator */

6      if $f(\pi )>f({\pi }^{\star })$ then

7          ${\pi }^{\star }\leftarrow \pi$, $NoImprove\leftarrow 0$

8       end

9     else if $vonlyconflictswithu\in \pi andhasnotbeeninvolvedinanySWAPoperation$ then

10      Remove u from π and insert v into π     /* SWAP operator */

11   else

12     Calculate NM(v) according to Eq. (7)

13     if $NM(v)\ne \mathrm{\varnothing }andvhasnotbeeninvolvedinanySWAPoperation$ then

14      Randomly select a vertex u from NM(v)

15      Remove v from π and insert u into π      /* SWAP operator */

16     else if $f(\pi )-1>f({\pi }^{\star })-\delta$ then

17      Remove v from π     /* DROP operation */

18   end

19 end

20 return $\pi ,{\pi }^{\star },NoImprove$

DOI: 10.7717/peerj-cs.1245/table-13

(1) If $v$ is outside $\pi$ (i.e., $v\in C_0$), the candidate move operator set consists of INSERT and SWAP. INSERT is applied first as it improves the solution quality. Then SWAP is applied, which keeps the solution quality unchanged. If neither operator can be applied, DTS just skips $v$. INSERT can be applied if the number of conflicts of $v$ in $\pi$ is 0 (i.e., $g(\pi \oplus <v,i_{-}(v),i_{+}(v)>)=0$). SWAP is applied if $v$ meets two conditions simultaneously: (1) $v$ is not involved in any SWAP operation already taken at the current round; (2) $v$ conflicts with just one vertex $u$ in $\pi$ (i.e., $g(\pi \oplus <v,i_{-}(v),i_{+}(v)>)=1$, and can only be swapped with $u$).

Similar to the INSERT operation in the greedy initialization, we adopt $g(\pi \oplus <v,i_{-}(v),i_{+}(v)>)$ for quick computation during the DTS stage. That is, we only have to update the number of conflicts of the operated vertex $v$ and its neighbor vertices after each move operation, i.e., the updated vertex set is $\{v\}\cup \{u:(u,v)\in Eandu\in C_0\}$. For the INSERT operation, we update the number of conflicts of the vertices neighboring to $v$ and not in $\pi$; for the SWAP operation, we update the number of conflicts of $v$, $u$, and all vertices neighboring to $u$ and $v$ not in $\pi$. Thus, the time complexity of INSERT and SWAP is $O(d_{max}^2)$, where $d_{max}$ denotes the largest degree of a vertex in the graph.

(2) If $v$ belongs to $\pi$, DROP and SWAP are the two candidate operators. SWAP is applied before DROP as SWAP does not degrade the solution quality while DROP does. If neither operation can be applied, the algorithm just skips $v$. SWAP can be applied only if $v$ satisfies two conditions simultaneously: (1) $v$ was not involved in any SWAP operation already taken at the current round; (2) the set $NM(v)\subset C_0$ of $v$ is non-empty, which is defined as

(7) $$NM(v)=\{u:g(\pi \setminus \{v\}\oplus <u,i_{-}(u),i_{+}(u)>)=0,c(u,v)=1\}.$$

The vertex $u$ that is to be swapped with $v$ is a random vertex in $NM(v)$.

DROP can be applied if it makes the number of the vertices in $\pi$ still above the threshold determined by $f({\pi }^{\star })-\delta$ after the DROP operation, where ${\pi }^{\star }$ is the best recorded solution and $\delta$ (a small positive integer) is a parameter. For the DROP operation, we need to update the number of conflicts of $v$ and all vertices neighboring to $v$ and not in the solution $\pi$. The time complexity of DROP is $O(d_{max}^2)$.

Figure 6 shows an example of the dynamic thresholding search stage. To explain this stage, we consider the solution in Fig. 5B as the input solution. For Fig. 5B, $C_0=\{a,b,d\}$, $U_0=\pi =\{e,c\}$ and $V_0=\{a,b,c,d,e\}$. Suppose that the vertices in $V_0$ are randomly shuffled into $\{a,e,d,b,c\}$. As shown in Fig. 6A, since the first vertex $a$ is outside $\pi$, INSERT and SWAP are the two candidate operators. INSERT is chosen to be applied before SWAP. The INSERT operator cannot be used since the number of conflicts of $v$ is not 0. However SWAP can be applied since $c$ only conflicts with $a$ and $a$ is not in the tabu list. As shown in Fig. 6B, for the second vertex $e\in \pi$, SWAP and DROP are the two candidate operators to be considered. SWAP is applied before DROP. SWAP can be applied since in this case $NM(e)=\{b,d\}$ and $v$ is not forbidden by the tabu list. Thus the second vertex $e$ is swapped with a random vertex in $NM(e)$, such as $d$. After that, the remaining vertices in $V_0$ are evaluated in the same way, while no operators can be applied to them. As a result, the improved solution is $\pi =\{d,a\}$.

Algorithm 5: :

The descent search

Input: Reduced graph $G=(V_0,E_0)$, solution π, best solution ${\pi }^{\star }$, the iterations without improvement NoImprove

Output: Solution π, best solution ${\pi }^{\star }$, the iterations without improvement NoImprove

1 while $N_I(\pi )\ne \mathrm{\varnothing }$ do

2   Randomly select a vertex u from $N_I(\pi )$ and insert it into π

3   if $f(\pi )>f({\pi }^{\star })$ then

4     ${\pi }^{\star }\leftarrow \pi$, $NoImprove\leftarrow 0$

5   end

6 end

7 return $\pi ,{\pi }^{\star },NoImprove$

DOI: 10.7717/peerj-cs.1245/table-14

The descent search stage

To complement the DTS stage where both equivalent and worse solutions are accepted, the descent search stage is subsequently applied to perform a more intensified examination of candidate solutions. Basically, this stage iteratively selects a conflict-free vertex and inserts it into the solution until such a vertex does not exist anymore.

Figure 7 shows an example of the descent search stage. To explain this stage, we consider the solution in Fig. 6B as the input solution, where $C_0=\{c,b,e\}$, $U_0=\pi =\{d,a\}$ and $V_0=\{a,b,c,d,e\}$. Through computation, $N_I(\pi )=\{e\}$. As shown in Fig. 7, the vertex $e$ from $N_I(\pi )$ is directly inserted into $\pi$ and the best solution ${\pi }^{\star }$ is updated to $\{d,e,a\}$.

Perturbation procedure

As described in “The Dynamic Thresholding Search Stage”, the threshold search accepts worse solutions that are within a certain quality threshold from the current solution, which relieves the search from the local optimum trap. However, there is a possibility that this strategy may fail. Therefore, we introduce a perturbation strategy that comes into effect when the search falls into a deep stagnation (i.e., the best solution does not change after $\omega$ consecutive local search runs). The perturbation strategy incorporates a learning mechanism that gathers move frequencies information from the local search, which is then advantageously used to guide the perturbation.

Algorithm 6 displays the perturbation procedure, which is decomposed into two steps:

Step 1: Choose and sort L vertices in $\pi$. Choose L vertices in $\pi$ with the highest move frequencies and sort them in a non-increasing order of the frequencies (lines 1–3, Algorithm 6). The move frequency of each vertex $v$ is the number of times that $v$ has been moved during the local search, which is initially set to 0, and increases by 1 each time $v$ is moved from one set to another.

Step 2: Drop and insert the to-be-perturbed vertices. Each vertex $v$ $(v\in A)$ is dropped, whose order $j$ in A is recorded (line 6, Algorithm 6). If $j>\lceil {\beta }_1\times (|\pi |+1)\rceil$ and $N_I(\pi )\ne \mathrm{\varnothing }$, randomly select a vertex $u$ from $N_I(\pi )$ and insert it into $\pi$ (lines 8–9, Algorithm 6). Recall that $N_I(\pi )$ represents the set of vertices in $C_0$ satisfying the condition $g(\pi \oplus ⟨v,i_{-}(v)>,i_{+}(v)⟩)=0$.

Figure 8 shows an example of the learning-based perturbation applied to a local optimal solution as shown in Fig. 7, where $C_0=\{c,b\}$ and $U_0=\pi =\{d,e,a\}$. Suppose $L=2$ and the chosen vertices is sorted as $A=\{e,a\}$ according to the move frequencies. The first vertex $e$ is dropped, which leads to an intermediate perturbed solution $\pi =\{d,a\}$. Then, the second vertex $a$ is also dropped. Since the order of $a$ is 2, which is more than $\lceil {\beta }_1\times (|\pi |+1)\rceil =1$ vertex, and $N_I(\pi )\ne \mathrm{\varnothing }$, the vertex $b$ is selected randomly from $N_I(\pi )$ and inserted into $\pi$, giving the perturbed solution $\pi =\{b,d\}$.

Recovery procedure

This is a reversed procedure of the reduction procedure. It restores the original graph $G=(V,E)$ from the reduced graph $G=(V_0,E_0)$ by adding back the removed vertices $U_r$ and the critical vertices $C_r$. Levy & Low (1988) indicates that the FVS of the original G is $C=C_0\cup C_r$. Figure 9 depicts an example that shows how the minimum FVS is determined. In the reduction procedure, $C_r=\{h\},U_r=\{f,g,i,j\}$. After the search stage, $V_0=\{a,b,c,d,e\}$ where $C_0=\{b,c\},U_0=\{a,d,e\}$. After the recovery procedure, the FVS of the original G is $C=C_0\cup C_r=\{b,c,h\}$.

Computational complexity and discussion

We consider first the greedy initialization procedure consisting of two stages. The first stage is to initialize the array $N_I(\pi )$, which can be realized in $O(|V_0|)$. The complexity of updating $N_I(\pi )$ is $O(d_{max}^2)$. The second stage is to construct the initial solution $\pi$, which is bounded by $O(|\pi |\times d_{max}^2)$, and $|\pi |$ is the size of $\pi$. Therefore, the time complexity of the greedy initialization procedure is $O(|V_0|+|\pi |\times d_{max}^2)$.

Next, the local search and perturbation procedures in the main loop of IDTS algorithm are considered. In each iteration of the local search, the dynamic threshold search and the descent search stages are performed alternately. The former is realized in $O(|\pi |\times d_{max}+|V_0|)$, and the latter in $O(|V_0|)$. Thus, the complexity of the local search procedure is $O(K_1\times (|\pi |\times d_{max}+|V_0|))$, where $K_1$ is the number of iterations of the local search. Then, the perturbation procedure can be achieved in $O(|\pi |\times ({\beta }_1+{\beta }_2\times d_{max}^2))$, which is much smaller than that of the local search. Therefore, the complexity of one iteration of the main loop of IDTS algorithm is $O(K_1\times (|\pi |\times d_{max}+|V_0|))$, and that of SA is $O(K_2\times (d_{max}^2+|V_0|))$, where $K_2$ is the number of iterations during each temperature period. Therefore, it can be seen that the two complexities are of the same order of magnitude.

Benchmark instances

We use 101 benchmark instances, which are classified into five categories. No optimal solutions are known for the instances of the first to forth categories, while optimal solutions are known for the instances of the fifth category.

• 1. The first category consists of 40 instances that are randomly generated by Pardalos, Qian & Resende (1998) using the FORTRAN random graph generator mkdigraph.f (http://mauricio.resende.info/data/index.html). The name of these instances is in the form of P|V | − |E|*, where $|V|\in \{50,100,500,1000\}$ is the number of vertices in the graph, and $|E|\in [100,30000]$ is the number of edges. Given the number of vertices and edges, a graph is built by randomly selecting |E| pairs of vertices as two endpoints of a directed edge. These instances are largely tested in the literature on the FVSP (Galinier, Lemamou & Bouzidi, 2013; Zhou, 2016).

• 2. The second category is composed of 10 random directed graphs, which are generated in the same way as the first category while the in-degree and out-degree of each vertex are no more than 10. These instances have R|V | − |E|* in their names. The number of vertices |V | is in the interval [100, 3,000], and the number of edges |E| in [500, 15,000].

• 3. The third category contains 10 artificially generated scale-free instances. These instances are generated by this work through the “powerlaw_cluster_graph” function of the “NetworkX” package, which is based on the algorithm proposed by Holme & Kim (2002). These instances are named as S|V | − |E|*, where |V | is in the interval [500, 3,000] and |E| is in [4,900, 29,900].

• 4. The fourth category is composed of 10 real-world instances from the Stanford large network dataset collection (http://snap.stanford.edu/data/). Nine of these instances are snapshots of the Gnutella peer-to-peer file sharing network. The remaining instance is a temporal network representing Wikipedia users editing each other’s Talk page. The number of vertices |V | is in the interval [6,301, 1,140,149], and the number of edges |E| in [20,777, 7,833,140].

• 5. The fifth category is composed of the 31 classical (easy) ISCAS89 benchmark instances which are from digital sequential circuits (Brglez, Bryan & Kozminski, 1989). These instances have s* in their names, where the number of vertices is in the range of [3, 1,728], and the number of edges in the range of [4, 32,774]. These instances, whose optima are known, are largely tested in the literature on the FVSP (Levy & Low, 1988; Lin & Jou, 1999; Orenstein, Kohavi & Pomeranz, 1995).

Experiment settings

The IDTS algorithm is programmed in C++ and compiled by GNU g++ 4.1.2 with the -O3 flag. Experiments are carried out on a computer with an Intel(R) Core(TM)2 Duo CPU T7700 2.4 GHz processor with 2 GB RAM running Ubuntu CentOS Linux release 7.9.2009 (Core).

Comparison with state-of-the-art results

Comparison of the results on the second-category instances

Table 4 shows the comparative results between IDTS and the reference algorithms on the 10 instances of the second-category. In terms of the best results, IDTS dominates Re-Red+SA by obtaining better values for all instances, and BPD by obtaining 7 better, and 3 equal results. On the other hand, IDTS significantly outperforms Re-Red+SA and BPD in terms of the worst and average results by obtaining better or equal results for all instances (except for R3000-15000). The small p-values (<0.05) indicate that there are significant differences between our best results and those of the two reference algorithms Re-Red+SA (p-value = 1.60E−3) and BPD (p-value = 8.20E−03).

Table 4:

Comparative results of IDTS with state-of-the-art algorithms on the 10 benchmark instances of the second category in normal test (30 independent runs).

Instance	|V |	|E|	Re-Red+SA				BPD (Zhou, 2016)				IDTS				${\mathrm{\Delta }}_1$	${\mathrm{\Delta }}_2$
			Best	Worst	Avg	$t(s)$	Best	Worst	Avg	$t(s)$	Best	Worst	Avg	$t(s)$
R100-500	100	500	34	41	36.6	0.00	30	33	31.0	0.32	30	30	30.0	0.03	−4	0
R200-1000	200	1,000	66	74	70.3	0.01	58	62	60.2	2.04	57	57	57.0	0.37	−9	−1
R500-2500	500	2,500	143	167	156.8	3.09	139	144	141.1	13.86	136	138	136.7	15.67	−7	−3
R800-4000	800	4,000	264	287	274.8	0.47	220	227	224.1	35.78	216	218	217.1	37.9	−48	−4
R1000-5000	1,000	5,000	302	332	319.0	11.64	262	276	267.6	59.05	262	267	264.7	30.14	−40	0
R1250-7500	1,250	7,500	474	518	500.5	40.74	406	411	409.1	139.70	403	411	406.6	98.62	−71	−3
R1500-9000	1,500	9,000	578	625	599.6	68.17	480	490	484.3	200.43	477	486	481.7	169.66	−101	−3
R1750-10500	1,750	10,500	658	727	694.9	131.02	561	570	566.0	270.63	561	569	564.2	287.72	−97	0
R2000-12000	2,000	12,000	726	808	763.6	163.09	640	650	644.5	351.42	638	648	643.5	340.99	−88	−2
R3000-15000	3,000	15,000	835	874	850.9	508.32	786	801	791.5	326.78	785	801	792.7	567.09	−50	−1
#Better			0	0	0		0	0	1
#Equal			0	0	0		3	2	0
#Worse			10	10	10		7	8	9
p-value			1.60E−03	1.60E−03	1.60E−03		8.20E−03	4.70E−03	1.14E−02

DOI: 10.7717/peerj-cs.1245/table-4

Note: The bold numbers in the table highlight the dominating results between the compared algorithms in terms of Best, Worst and Avg values.

Furthermore, Fig. 10 summarizes the performance of the IDTS algorithm with that of the Re-Red+SA and BPD algorithms on these instances. Figure 10A presents the relationship between the number of vertices and the best FVS size (the best objective value over 30 runs). Figure 10B shows the relationship between the number of vertices and the average computation time. One observes that the FVS size increases linearly while the average computation time increases exponentially with the increase of the number of vertices.