When does the Physarum Solver Distinguish the Shortest Path from other Paths: the Transition Point and its Applications
Yusheng Huang, Dong Chu, Joel Weijia Lai, Yong Deng, Kang Hao Cheong
IIEEE TRANSACTIONS 1
When does the Physarum Solver Distinguish theShortest Path from other Paths: the Transition Pointand its Applications
Yusheng Huang, Dong Chu, Joel Weijia Lai, Yong Deng, Kang Hao Cheong
Abstract —Physarum solver, also called the physarum poly-cephalum inspired algorithm (PPA), is a newly developed bio-inspired algorithm that has an inherent ability to find the shortestpath in a given graph. Recent research has proposed methodsto develop this algorithm further by accelerating the originalPPA (OPPA)’s path-finding process. However, when does thePPA ascertain that the shortest path has been found? Is therea point after which the PPA could distinguish the shortest pathfrom other paths? By innovatively proposing the concept of thedominant path (D-Path), the exact moment, named the transitionpoint (T-Point), when the PPA finds the shortest path can beidentified. Based on the D-Path and T-Point, a newly acceleratedPPA named OPPA-D using the proposed termination criterionis developed which is superior to all other baseline algorithmsaccording to the experiments conducted in this paper. The validityand the superiority of the proposed termination criterion is alsodemonstrated. Furthermore, an evaluation method is proposed toprovide new insights for the comparison of different acceleratedOPPAs. The breakthrough of this paper lies in using D-pathand T-point to terminate the OPPA. The novel terminationcriterion reveals the actual performance of this OPPA. ThisOPPA is the fastest algorithm, outperforming some so-calledaccelerated OPPAs. Furthermore, we explain why some existingworks inappropriately claim to be accelerated algorithms is infact a product of inappropriate termination criterion, thus givingrise to the illusion that the method is accelerated.
Index Terms —Physarum polycephalum inspired algorithm,Bio-inspired algorithm, Physarum Solver, Shortest path problem
I. I
NTRODUCTION
Bio-inspired algorithm (BA) has attracted attention andintrigued many due to its wide range of applications for along time. Numerous kinds of BAs, such as the differentialevolution algorithm [1], grey wolf optimizer [2], artificialbee colony algorithm [3], genetic algorithm [4], ant colonyalgorithm (ACO) [5], and whale optimization algorithm [6]have been developed. BAs have played important roles in solv-ing optimization problems, such as large-scale multiobjectivedetection and optimization [7], [8], parameter estimation [9],global numerical optimization [10], feature selection [11], nu-merical function optimization [3], and boosting support vectormachine [12], just to name a few. In the past decade, a new BA,
Y. Huang, D. Chu and Y. Deng are with the Institute of Fundamen-tal and Frontier Science, University of Electronic Science and Technol-ogy of China, Chengdu, 610054, China (e-mail: [email protected],[email protected]).J.W. Lai and K.H. Cheong is with Science and Math Cluster, SingaporeUniversity of Technology and Design (SUTD), S487372, Singapore (e-mail:kanghao [email protected]). named the physarum polycephalum inspired algorithm (PPA)[13], [14], also called the physarum solver, has attracted greatattention.The physarum polycephalum, a kind of slime mold, is asingle-celled amoeboid organism [13]. In 2000, Nakagaki etal. [15] observed that the amoeboid organism would gradu-ally change its shape to improve its efficiency of foraging,which leads to the result that its body would only cover theshortest path in the mazes after some time. The maze-solvingability of the slime mold was then considered as a kind ofprimitive intelligence [16]. It is believed that the maze-solvingmechanism which is also the path-finding mechanism wasrelated to the contraction waves in the organism [17]. Thecontraction waves influence the thickness of the tube, resultingin a positive mechanism. If the flow through a certain tubepersists or increases for a certain period, the tube will becomethicker; otherwise, the tube will contract [18].In 2007, a mathematical model of the physarum poly-cephalum, which is the original physarum polycephalum in-spired algorithm (OPPA) was proposed [14]. In the OPPA, thetube-liked body of the physarum polycephalum is consideredas a network with the growing point and the food sourceserving as the starting point and ending point, respectively;the Poiseuille flow is adopted to model the flow flowingthrough the tube; the conservation law of flux (also calledcontinuity of flux) is obeyed using network Poisson equation[13]. At the core of the OPPA is the use of a novel adaptiveequation to model the dynamics of tube thickness. By using thetechniques mentioned above, the OPPA perfectly models thepath-finding behavior of the amoeboid organism and was thenapplied to solve the minimum-risk path-finding problems [14]and adaptive network design problems [13]. The effectivenessof the OPPA has been mathematically proven. Bonifaci etal. provided a mathematical proof of the convergence of theOPPA on the shortest path problems [19] and provided thetime complexity bounds [20]. Karrenbauer et al. developedon the former proof and proposed a more general proof forthe effectiveness of the OPPA [21].Many relevant research have been proposed after the initialintroduction of the OPPA. Some studies adopt the PPA toimprove the performance of other BAs. For example, PPA isutilized to initialize the ant colony algorithm to solve someNP-Hard problems [22]; Gao et al. combined some traditionalBAs with PPA in network community detection; a combinedPPA framework with the genetic algorithm is proposed tohandle the traveling salesman problem [23], just to name a r X i v : . [ c s . N E ] J a n EEE TRANSACTIONS 2 a few. Other studies aim to expand the PPA’s applications.For example, Sun et al. [24] proposed a physarum-inspiredapproach to identify sub-networks for drug re-positioning; Xu et al. [25] developed a multi-sink-multi-source PPA to tacklethe traffic assignment problem; Tsompanas et al. [26], andZhang et al. [27] also worked on finding equilibrium state oftraffic assignments using PPA; Jiang et al. [28] applied thePPA in routing protocol design; Song et al. [29] designed aPPA-based optimizer for minimum exposure problem, etc.Besides the research mentioned above, some researchers arealso interested in improving the convergence rate of the OPPA.The computational complexity for solving the network Poissonequation using PPA is O ( n ) . While it is polynomial time, thetypical number of variables is of a much higher order, thusresulting in an inefficient allocation of time resource beforeconverging to the shortest path [18]. Thus, several studiesrelated to accelerated OPPA have been developed. Zhang et al. [18] proposed an improved OPPA by introducing theconcept of energy. In the improved PPA, it is assumed thatthe foraging of the physarum polycephalum consumes energyand the law of the conservation of energy should be satisfied.In another work, the OPPA is accelerated by eliminating somenear vanished tubes of the OPPA [30]. Wang et al. [31]believed that the tubes in the shortest path would be morecompetitive than others when OPPA converges. Thus theyproposed an anticipation mechanism to terminate the OPPAearlier. Cai et al. [32] combined the OPPA with the Bayesianrule to achieve a higher convergence rate and also proposeda method to tackle negative-weighted edges in the OPPA.Gao et al. [33] developed an accelerated OPPA by excludingthe inactive nodes and collecting the near-optimal paths. Theabove researchers share the same motivation—while the OPPAis effective, its efficiency can be greatly improved; henceaccelerating the OPPA is necessary. These work have spedup the OPPA using different methods. Yet, despite significantachievements being made by previous studies, there remainsome questions that need to be answered:1) In the case of ensuring accuracy, is there an inherentlimit so that no matter how the termination criteriachanges, the OPPA cannot increase the convergencerate?
Both [31] and [33] stated the necessity of defining anappropriate stopping condition for the OPPA and claimedthat the proposed solutions could terminate the OPPAearlier. This leads to a secondary question, how muchearlier? Is there a limit for the early termination, prior towhich the OPPA have not converged to the shortest path?2)
Is there a better evaluation method that could moreintuitively evaluate how much improvement the accel-erated methods are making?
All of the studies mentionedabove used either the running time or the number ofiterations to compare different accelerated OPPAs. How-ever, the two metrics are highly dependent on the chosenstopping criteria, in other words, if the evaluation methodis dependent on the stopping criteria, the comparison maynot be objective. Hence, is there an evaluation method thatis not dependent to the stopping criteria of the OPPA?The two questions above ultimately lead to the main question: when does the OPPA distinguish the shortest path fromthe other paths?
The contributions of this work are listed below: • We have defined two concepts, i.e., the dominant path andthe transition point of the OPPA. By combining the useof these two concepts, we identify the exact moment thatthe OPPA starts to detect the shortest path for the firsttime. These concepts answer the first and main questions. • Combining the defined concepts, we propose a new stop-ping criterion for the OPPA. By combining the proposedcriterion and the OPPA, a new algorithm named theOPPA-D is developed. According to the experimentalresults in Section IV-A, the OPPA is the fastest among allof the tested algorithms. The validity and the superiorityof the proposed criterion are also demonstrated in SectionIV-B. • By adopting the concept of the transition point, wehave developed a novel method to compare differentaccelerated OPPAs more objectively. By doing so, thesecond question can be answered. Using the proposedmethods, we compare our algorithm with existing acceler-ated OPPAs in Section IV-C where we also present somekey findings.Broadly, the paper is developed as follows: Section II intro-duces the background and development of the PPA; Section IIIprovides the defined concepts and the proposed algorithm; theproposed algorithm is then experimentally tested in SectionIV; Section IV also provides the proposed evaluation methodand the corresponding analysis. Section V discusses andconcludes the paper.II. B
RIEF INTRODUCTION OF THE ORIGINAL
PPAThis section is mainly based on [13], [14], where moredetails could be found. The pseudo-code of the OPPA isprovided in Section II-B for the readers’ reference.
A. The mathematical model of the physarum polycephalum
As mentioned, the slime mold’s body would form a tube-like network when foraging. Consider the network as a graph G ( N, E ) where N is the set of nodes and E represents theset of edges. Note that G is assumed to be an undirectedconnected graph without cycles and negative-weight edges,except where specifically mentioned. Other abstract mappingsare as follows: the tubes of the network are treated as edges in G ( N, E ) ; the junctions between the tube are regarded as nodesin G ( N, E ) ; the growing point of the slime mold is mappedto the ending node of the G ( N, E ) while the food source isthe starting node. For simplicity, in this paper, the scenariowith one starting node and one ending node is considered.For scenarios with multiple ending nodes please refer to [34]and as for scenarios with multiple starting nodes please referto [25].In this paper, N and N denotes starting and the end-ing nodes, respectively; the other nodes are labeled as N , N , N , · · · ; E ij represents the edge between nodes N i and N j ; Q ij represents the flux flowing through the edge E ij from node N i to node N j . EEE TRANSACTIONS 3
In [14], the flow in the slime mold’s network is regardedas laminar. Thus, the Hagen-Poiseuille equation is observed,giving us Q ij = D ij L ij ( p i − p j ) , (1)where p i denotes the pressure at node N i ; D ij represents theconductivity of the edge E ij and L ij is the length of the edge.The whole network is driven by the inflow at the food source(the starting node). Thus, at N , the source term of flux is (cid:88) i (cid:54) =1 Q i = IN (2)where IN is the inflow of the network.The outflow of the network is considered equal to the inflowwhich is IN . Then, at the ending node ( N ) of the graph,we have (cid:88) i (cid:54) =2 Q i = − IN (3)The inflow and outflow at other nodes should be balanced,that is (cid:88) j (cid:54) =1 , i (cid:54) = j Q ij = 0 . (4)By combining (1) to (4), the network Poisson equation ofthe graph G ( N, E ) is (cid:88) i (cid:54) = j D ij L ij ( p i − p j ) = + IN for j = 1 , − IN for j = 2 , otherwise . (5)All p i can be calculated by solving the network Poissonequation, i.e., (5), by setting p to . After obtaining p i , Q ij can be determined using (1). However, we still need D ij tosolve (5). We now discuss the method to update D ij at eachiteration which is the main feature of the OPPA.In [13], [14], the following adaptation equation is adoptedto model the dynamics of the thickness of tubes: ddt D ij = f ( | Q ij | ) − αD ij , (6)where f ( | Q ij | ) = | Q ij | and α = 1 are typically used. Furtherdiscussion on different parameter settings of the adaptationequation can be found in previous literature [14], [21]. Thisadaptation equation suggests that the conductivity of tubetends to decrease exponentially while it increases linearly withflux along this tube.In order to implement the adaptation equation, we discretize(6) by performing linear approximations to get D n +1 ij − D nij ∆ t = | Q nij | − D n +1 ij , (7)where Q nij is the flow flowing through edge E ij at n thiteration; ∆ t is normally considered as (please refer to [33]for various other considerations); D nij and D n +1 ij representsthe conductivity of edge E ij at n -th and ( n + 1) -th iteration,respectively. A more concise form of the above formula is D n +1 ij = | Q nij | + D nij . (8) B. The pseudo-code of the OPPA
Having introduced the mathematical model for thephysarum solver, we will now provide the pseudo-code ofthe OPPA. For the rest of the paper, ‘OPPA’ will refer tothe following algorithm which uses the original PPA to solvethe shortest path problem. The pseudo-code of the OPPA isdemonstrated in Alg.1.
Algorithm 1
The OPPA [13]
1: //Initialization part2:
Input: the statistics of graph G ( N, E ) and the corresponding matrix ofthe weight of each edges W .3: D ij ⇐ . ∀ E ij ∈ E, ) // Initialize the conductivity of each edge.4: Q ij ⇐ ∀ E ij ∈ E ) // Initialize the flow of each edge.5: L ij ⇐ W ij ( ∀ E ij ∈ E ) //The length of the edges equal to their weights.6: p i ⇐ ∀ i = 1 , , · · · , N ) // Initialize the pressure at each node.7: count ⇐ // Initialize the counting variable.8: //Iterative part9: repeat p ⇐ // The pressure of the starting node (Node ) is set to 0.11: Calculate the pressure of all the nodes using (5): (cid:88) i ∈ N D ij L ij ( p i − p j ) = + IN for j = 1 , − IN for j = 2 , otherwise.
12: Calculate the flux using (1): Q ij ⇐ D ij · ( p i − p j ) /L ij .
13: Calculate the conductivity of the iteration using (8): D n +1 ij = | Q nij | + D nij . count ⇐ count + 1 until The given termination criterion is met.16:
Output:
The flow matrix Q . Given a graph, after the execution of the OPPA, accordingto the flow matrix Q , one would find a path that contains themost amount of inflow in the graph, this is the shortest path.III. D EFINING CONCEPTS AND THE PROPOSED
OPPA-D
ALGORITHM
A. The dominant path
Dominant path (D-Path):
Given a flow matrix generated bythe PPA-based algorithms, the D-path is the path found by thefollowing iterative procedures:i. Set up a node list
List
Node initialized as an empty list.The starting node N of the graph would be the firstelement of List
Node .ii. To find the next element of
List
Node in the graph, find allthe nodes and their corresponding edges that are connectto the current last element of
List
Node .iii. Among the edges found in the second step, find the edgethat contains the maximum flow. Since one of the twonodes that relate to this edge is the last element of thenode set
List
Node , set the other node as the next elementof
List
Node .iv. Delete the current last element of
List
Node and its corre-sponding edges from the graph. Push the next element of
List
Node into the list
List
Node . This element would bethe last element of
List
Node in the next iteration.
EEE TRANSACTIONS 4 v. Repeat Step ii to Step iv until the next element we foundin Step iv is the ending node of the graph. The nodes in
List
Node form one of the paths in the graph. This pathis defined as the
D-Path .Given a flow matrix generated by the PPA-based algorithms,the D-Path can be generated following the above procedures,and the length L D − P ath of the path can be calculated.
B. The transition point of the physarum solver
With the definition of D-Path, it is interesting to note thatafter the convergence of the OPPA, the algorithm will directmost (or even all) of the inflow to the shortest path; thus,when the OPPA converges, the D-Path generated is exactlythe shortest path of the graph, and the L D − P ath will be equalto the length of the shortest path. One can also deduce thatthe process of the OPPA finding the shortest path is indeedthe process of OPPA gradually converging the D-Path to theshortest path.However, the means to quantify the convergence of theOPPA is still an open problem. The convergence of the OPPAis typically regarded as the iteration when the conductivity ofeach edge in the graph does not change with significantly withthe subsequent iterations. Thus, the stopping (or convergenceor termination) criterion is usually set to (cid:80) i (cid:80) j (cid:54) = i D ij ≤ (cid:15) . In[33] (cid:15) is set to − while in [18], [30] it is set to − . Thereis still no conclusive method or consensus for determining theparameter (cid:15) .In most of the current studies, we allow the algorithm torun until it achieves the pre-determined stopping criterion.This would often result in long computational time; yet itremains a widely implemented method for deciding whetherthe OPPA finds the shortest path. The question remains, doesthe OPPA need such a long computational time to achieveconvergence? With our definition of D-Path, the answer is no . The OPPAdistinguishes the shortest path before the stopping criterion ismet. To demonstrate this, four complete graphs, i.e., ’Instance1’ to ’Instance 4’, are randomly generated with five nodesand the length of edges varying from 1 to 10000. The OPPAis then implemented in the four instances to find the shortestpath. The termination criterion is set to (cid:80) i (cid:80) j (cid:54) = i D ij ≤ − .At each iteration, the L D − P ath is recorded.Fig.1 demonstrates the results. According to Fig.1, in allof the instances, the L D − P ath fluctuates at the first fewiterations, and then decreases rapidly, after which it convergesto the length of the shortest path. Of importance is that theconvergence of L D − P ath happens long before the OPPA’stermination. Thus, from Fig.1, one could tell that the OPPAdistinguishes the shortest path long before the terminationcondition is met.
Transition point (T-Point):
The T-point of the PPA-basedalgorithms is defined as the iteration when the L D − P ath startsto converge to the length of the shortest path. For example, thefour points with green arrows pointing at in Fig.1 are the T-Points of the OPPA in each of the four instances, respectively.Given the definition of T-Point, the main question proposedin Section I, i.e., “When does the OPPA distinguish the
Iteration
Leng t h o f t he do m i nan t pa t h Instance1Instance2Instance3Instance4
Fig. 1. The transition point. shortest path from the other paths?”, can be answered. Wewould like to provide our opinion: after the transition point,the PPA-based algorithms start to distinguish the shortest pathfrom the others. For the first question in Section I, i.e., “Inthe case of ensuring accuracy, is there an inherent limit, sothat no matter how the termination criteria changes, the OPPAcannot increase the convergence rate?”, again, we would liketo propose a possible answer: the inherent limit would be theT-Point, before which the OPPA does not achieve the shortestpath.
C. The proposed algorithm: OPPA-D
Further developing from the previous sub-sections, we nowintroduce our application of the proposed concepts, i.e., theOPPA-D algorithm.The OPPA-D algorithm builds on the OPPA with a newlyproposed termination criterion. This criterion is defined asfollows: Execute the OPPA until the L D − P ath does not changefor K iteration. K is a pre-defined parameter. Thus, after the L D − P ath remains steady for K iterations, the occurrence ofthe T-Point is assumed, which also implies the convergence ofthe OPPA.The pseudo-code of the OPPA-D algorithm is providedin Alg.2. For purpose of illustration, K is set to . Thealgorithm’s sensitivity to K will be subsequently tested inSection IV-B. Algorithm 2
The OPPA-D
1: //Initialization part2:
Input: the statistics of graph G ( N, E ) , the corresponding matrix of theweight of each edges W , and the pre-defined parameter K .3: D ij ⇐ . ∀ E ij ∈ E, ) // Initialize the conductivity of each edge.4: Q ij ⇐ ∀ E ij ∈ E ) // Initialize the flow of each edge.5: L ij ⇐ W ij ( ∀ E ij ∈ E ) //The length of the edges equal to their weights.6: p i ⇐ ∀ i = 1 , , · · · , N ) // Initialize the pressure at each node.7: count ⇐ // Initialize the counting variable.8: Length D − Path ⇐ // Variable used to record the length of D-Path.9: COUNT ⇐ // Variable used to count the number of iterations thatthe Length D − Path remains unchanged.10: //Iterative part11: while true do
EEE TRANSACTIONS 5 p ⇐ // The pressure of the starting node (Node ) is set to 0.13: Calculate the pressure of all the nodes using (5): (cid:88) i ∈ N D ij L ij ( p i − p j ) = + IN for j = 1 , − IN for j = 2 , otherwise.
14: Calculate the flux using (1): Q ij ⇐ D ij · ( p i − p j ) /L ij .
15: Calculate the conductivity of the iteration using (8): D n +1 ij = | Q nij | + D nij .
16: Find the D-Path and calculate its length
Length countD − Path accordingto Section III-A.17: if Length countD − Path == Length count − D − Path then
COUNT = COUNT + 1 if COUNT ≥ K then break end if else COUNT = 0 end if count ⇐ count + 1 end while Output:
The shortest path length
Length countD − Path . IV. E
XPERIMENTS & A
NALYSIS
In this section, we present our experimental methodologiesand our findings. In particular, we compare the proposedOPPA-D with other state-of-the-art accelerated OPPAs inrandomly generated complete networks, randomly generatedsmall-world networks, and real-world networks, to demon-strate its effectiveness and efficiency in Section IV-A. Since theOPPA-D is a modified algorithm of the OPPA with the newlyproposed termination criterion, we will verify the validation ofthe proposed stopping criterion. Hence, the OPPA-D is furthercompared with other commonly adopted criteria in SectionIV-B. In Section IV-C, three accelerated OPPAs will be evalu-ated using a new transition-point based evaluation method. Allthe experiences are done through computer simulations usingMatlab R2018a on an Intel Core i5-8500 CPU (3GHz) with 8GB RAM under Windows 10.
A. Comparing the OPPA-D with other accelerated OPPAs
Data Set 1:
There are three data sets. 1) ‘Da-Com-1’ to ‘Da-Com-8’: Eight complete graphs are randomly generated withsize ranging from 100 to 5000, and length of the edges varyingfrom 1 to 10000. 2) ‘Da-SW-1’ to ‘Da-SW-8’: Eight small-world networks are randomly generated with size rangingfrom 100 to 5000, and length of the edges varying from 1to 10000. All the small-world networks are created using aMatlab function ’WattsStrogatz’ . 3) ‘Da-R-1’ and ‘Da-R-2’:Two real-world networks, i.e., the Sioux-Falls network andthe Anaheim network , are adopted. Some basic informationabout the data sets is provided in Table I. More details of the function ’WattsStrogatz’ could refer tohttps://ww2.mathworks.cn/help/matlab/math/build-watts-strogatz-small-world-graph-model.html?lang=en. The meta-data of networks could refer tohttps://github.com/bstabler/TransportationNetworks. TABLE IB
ASIC INFORMATION OF D ATA S ET Instance Name | N | a | E | a Instance Name | N | a Mean node degree Instance Name | N | a | E | a Da-Com-1 50 1.22E+03 Da-SW-1 50 6 Da-R-1 b
24 76Da-Com-2 100 4.95E+03 Da-SW-2 100 12 Da-R-2 b
416 914Da-Com-3 250 3.11E+04 Da-SW-3 250 30Da-Com-4 500 1.25E+05 Da-SW-4 500 60Da-Com-5 750 2.81E+05 Da-SW-5 750 90Da-Com-6 1000 4.99E+05 Da-SW-6 1000 120Da-Com-7 2000 2.00E+06 Da-SW-7 2000 240Da-Com-8 5000 1.25E+07 Da-SW-8 5000 600 a N is the number of nodes and E is the number of edges. b Instances ‘Da-R-1’ and ‘Da-R-2’are the Sioux-Falls network and theAnaheim network, respectively.TABLE IIA
VERAGE RUNNING TIME OF TESTED ALGORITHMS IN D ATA S ET Instacnes Running time (Second)OPPA( − ) OPPA( − ) EHPA APS OPPA-D(Ours)Da-Com-1 0.0020 0.0049 0.0073 0.0023 Da-Com-2 0.0223 0.0557 0.0131
Da-Com-4 0.4839 1.0926 0.3771 0.2293
Da-Com-5 2.7072 6.1812 1.0617 0.6322
Da-Com-6 3.5256 8.2404 2.8248 1.3938
Da-Com-7 5.7222 9.4358 20.0010 16.0612
Da-Com-8 50.1900 89.3060 444.2410 277.3507
Da-SW-1 0.0015 0.0035 0.0010 0.0021
Da-SW-2 0.0086 0.0195
Da-SW-4 0.2700 0.5818 0.2176 0.2157
Da-SW-5 1.1505 2.7874 0.8128 0.5841
Da-SW-6 0.9320 1.4951 1.1180 1.1162
Da-SW-7 20.1512 53.4956 27.4835 15.3748
Da-SW-8 121.5098 277.3977 620.7715 279.5441
Da-R-1 0.0012 0.0021 0.0028 0.0039
Da-R-2 0.8085 1.8604 0.3038 0.2507 a Values presented in the table are the average value of 15 trials.
Baseline algorithms:
The OPPAs with termination crite-ria (cid:15) ≤ − and (cid:15) ≤ − , i.e., the OPPA( − ) andOPPA( − ), are utilized as the baseline algorithms. Two otherstate-of-the-art accelerated OPPAs, i.e., the EHPA [30] and theAPS [33] are also compared against . Each algorithm will beevaluated 15 times for each graph.All the evaluated algorithms successfully solved for theshortest path in all 15 instances. The running time of themare presented in Fig.2, 3, 4. Tables II and III summarize therunning time and the number of iterations, respectively. The codes of the EHPA and the APS are available athttps://github.com/caigaoub/PhysarumOptimization. We thank the authors forgiving access to their source codes.TABLE IIIA
VERAGE NUMBER OF ITERATIONS OF TESTED ALGORITHMS IN D ATA S ET Instacnes Number of iterationsOPPA( − ) OPPA( − ) EHPA APS OPPA-D(Ours)Da-Com-1 26 58 12 16 12Da-Com-2 57 134 21 16 19Da-Com-3 27 49 16 16 18Da-Com-4 58 139 23 17 25Da-Com-5 134 337 26 26 25Da-Com-6 96 238 33 20 22Da-Com-7 31 55 19 16 19Da-Com-8 41 75 26 17 23Da-SW-1 24 47 8 16 11Da-SW-2 27 50 14 16 11Da-SW-3 41 89 20 16 17Da-SW-4 34 73 16 16 11Da-SW-5 62 149 25 18 20Da-SW-6 27 43 16 16 17Da-SW-7 116 305 29 16 22Da-SW-8 105 233 38 36 30Da-R-1 27 53 9 16 11Da-R-2 133 330 37 44 24 a Values presented in the table are the average value of 15 trials.
EEE TRANSACTIONS 6
OPPA(1e-2)OPPA(1e-5) APS EHPA OPPA-D10 -3 -2 R unn i ng t i m e ( s e c ond ) (a) ‘Da-Com-1’ instance . OPPA(1e-2)OPPA(1e-5) APS EHPA OPPA-D10 -2 -1 R unn i ng t i m e ( s e c ond ) (b) ‘Da-Com-2’ instance . OPPA(1e-2)OPPA(1e-5) APS EHPA OPPA-D0.020.030.040.050.060.070.08 R unn i ng t i m e ( s e c ond ) (c) ‘Da-Com-3’ instance . OPPA(1e-2)OPPA(1e-5) APS EHPA OPPA-D0.20.30.40.50.60.70.80.911.1 R unn i ng t i m e ( s e c ond ) (d) ‘Da-Com-4’ instance . OPPA(1e-2)OPPA(1e-5) APS EHPA OPPA-D10 R unn i ng t i m e ( s e c ond ) (e) ‘Da-Com-5’ instance . OPPA(1e-2)OPPA(1e-5) APS EHPA OPPA-D10 R unn i ng t i m e ( s e c ond ) (f) ‘Da-Com-6’ instance . OPPA(1e-2)OPPA(1e-5) APS EHPA OPPA-D468101214161820 R unn i ng t i m e ( s e c ond ) (g) ‘Da-Com-7’ instance . OPPA(1e-2)OPPA(1e-5) APS EHPA OPPA-D10 R unn i ng t i m e ( s e c ond ) (h) ‘Da-Com-8’ instance .Fig. 2. Experimental results in complete graphs. OPPA(1e-2)OPPA(1e-5) APS EHPA OPPA-D10 -3 R unn i ng t i m e ( s e c ond ) (a) ‘Da-Sw-1’ instance . OPPA(1e-2)OPPA(1e-5) APS EHPA OPPA-D10 -2 R unn i ng t i m e ( s e c ond ) (b) ‘Da-Sw-2’ instance . OPPA(1e-2)OPPA(1e-5) APS EHPA OPPA-D10 -2 R unn i ng t i m e ( s e c ond ) (c) ‘Da-Sw-3’ instance . OPPA(1e-2)OPPA(1e-5) APS EHPA OPPA-D10 -1 R unn i ng t i m e ( s e c ond ) (d) ‘Da-Com-4’ instance . OPPA(1e-2)OPPA(1e-5) APS EHPA OPPA-D10 R unn i ng t i m e ( s e c ond ) (e) ‘Da-Sw-5’ instance . OPPA(1e-2)OPPA(1e-5) APS EHPA OPPA-D0.60.70.80.911.11.21.31.41.5 R unn i ng t i m e ( s e c ond ) (f) ‘Da-Sw-6’ instance . OPPA(1e-2)OPPA(1e-5) APS EHPA OPPA-D10 R unn i ng t i m e ( s e c ond ) (g) ‘Da-Sw-7’ instance . OPPA(1e-2)OPPA(1e-5) APS EHPA OPPA-D10 R unn i ng t i m e ( s e c ond ) (h) ‘Da-Sw-8’ instance .Fig. 3. Experimental results in small-world networks. EEE TRANSACTIONS 7
OPPA(1e-2)OPPA(1e-5) APS EHPA OPPA-D10 -3 -2 R unn i ng t i m e ( s e c ond ) (a) ‘Da-R-1’ instance . OPPA(1e-2)OPPA(1e-5) APS EHPA OPPA-D10 -1 R unn i ng t i m e ( s e c ond ) (b) ‘Da-R-2’ instance .Fig. 4. Experimental results in real-world networks. In general, the OPPA-D is superior to the baseline al-gorithms in terms of efficiency. In the complete graphs, asrecorded in Fig.2 and Table II, the OPPA-D outperforms allother baseline algorithms except for ‘Da-Com-2’ where theAPS is slightly faster than the OPPA-D. In the other instances,the OPPA-D is 56.40%, 0.37%, 12.48%, 25.69%, 44.61%,40.58%, and 44.91% faster than the next fastest algorithm,respectively. In the small-world graphs, according to Fig.3 andTable II, the OPPA-D is more efficient than all the baselinealgorithms except for ‘Da-Sw-2’ where it ranks fourth. In otherinstances, the OPPA-D is 26.21%, 39.19%, 59.06%, 37.14%,35.76%, 75.07%, and 71.34% faster than the second-rankedalgorithm, respectively. In the real-world network, accordingto Fig.4 and Table II, the OPPA-D is the fastest algorithm. Itis 63.09% and 38.84% more efficient than the algorithm insecond place, respectively.Though the accelerated OPPAs, i.e., the EHPA and the APS,are faster than the OPPA( − ) and OPPA( − ) in mostinstances, they do not out-perform the OPPA( (cid:15) = 10 − ) inthe 2000-nodes and 5000-node complete graphs and 5000-node small-world network. One interesting finding in thisexperiment is that, according to Tables II and III, while thenumber of iterations is small it does not necessarily translateto a shorter running time. For example, the OPPA( − ) takes41 iterations to converge while the APS only take 17 iterationsfor ‘Da-Com-8’. However, the OPPA( − ) is approximatelyfive times faster than the APS. Therefore, when comparingdifferent accelerated OPPAs, the running time, instead of thenumber of iterations, would be a more appropriate evaluationmetric. B. Comparing different termination criteria
In the above experiments, we observe that the runningtime of the OPPA is highly sensitive to the current stoppingcriteria, i.e., (cid:80) i (cid:80) j (cid:54) = i D ij ≤ (cid:15) . According to Table II, theOPPA( − ) is more efficient than the OPPA( − ); bothalgorithms solve for the shortest-path in all the instancessuccessfully. Recall that the OPPA-D is fundamentally theOPPA with a newly proposed termination criterion. Thus, tovalidate the superiority of the proposed criterion, i.e. executingthe OPPA until the L D − P ath reamins steady for K iteration,we conduct further experiments. The evaluated criteria:
1) The traditional termination cri-terion with (cid:15) = 10 − , − , · · · , − are tested. 2) The proposed termination criterion with K = 5 , , , , arealso tested. Data Set 2:
We test the termination criteria in randomlygenerated complete graphs. The size of the graphs are 10, 100,and 500. For each size, 50 graphs are randomly generated withthe length of the edge varying from 1 to 10000.Tables IV and V report the experimental results. From TableIV, different settings towards parameters (cid:15) and K may leadto different success rates. One interesting thing of note is thatwhen the graph size is 100 or 500, as (cid:15) decreases, despite thedecrease in the number of failed attempts, the algorithm couldnot reach the traditional termination criterion with a low (cid:15) incertain cases. However, the proposed stopping criterion doesnot have the above issue.According to Table V, one could conclude that, from theperpective of running time, the traditional termination criterionis more sensitive to the parameter (cid:15) , compared with theproposed stopping criterion’s sensitivity to K . For example,when the graph size is 500, when (cid:15) decreases from − to − , the running time increases from 0.32 seconds to 2.40seconds (by 650%); however, when K increases from to ,the running time increases from 0.10 seconds to 0.29 seconds(by 190%).In general, in terms of of success rate, the traditional termi-nation criterion with (cid:15) = 10 − and the proposed terminationcriterion with K = 30 would be good choices, both of whichachieve 100% success rate for three different graph sizes(each of the graph size has 50 cases). After comparing ourexperimental results of the running time for both criteria, onewould find that the proposed criterion with K = 30 would bethe better choice. C. Analyzing the improvements made by accelerated OPPAs
In the previous sub-section, we demonstrate that the pro-posed criterion is better than the traditional criterion. However,the running time of the OPPA for both stopping criteria is stillinevitably dependent on the initial parameters. It is logical todeduce that the above findings still holds when these criteriaare applied to other accelerated OPPAs. Thus, in this section,we want to compare the accelerated OPPAs with the OPPAby excluding the termination criteria.
1) The proposed evaluation method:
Using the definedconcepts of the D-Path and T-Point, we propose our evaluationmethod:
EEE TRANSACTIONS 8
TABLE IVE
XPERIMENTAL RESULTS IN D ATA S ET UCCESS RATE . Size The traditional criteria ( (cid:15) ) The proposed criteria ( K )( − ) ( − ) ( − ) ( − ) ( − ) 5 10 15 20 3010 96%(2 |
0) 100%(0 |
0) 100%(0 |
0) 100%(0 |
0) 100%(0 |
0) 96%(2 |
0) 96%(2 |
0) 100%(0 |
0) 100%(0 |
0) 100%(0 | |
0) 100%(0 |
0) 96%(0 |
2) 92%(0 |
4) 88%(0 |
6) 60%(20 |
0) 88%(6 |
0) 90%(5 |
0) 96%(2 |
0) 100%(0 | |
0) 100%(0 |
0) 98%(0 |
1) 96%(0 |
2) 96%(0 |
2) 64%(18 |
0) 92%(4 |
0) 98%(1 |
0) 100%(0 |
0) 100%(0 | a The values in the table have a format ‘Value1%(Value2 | Value3)’, where Value1 is the success rate out of 50 cases, Value2 represents the number of thecases that the algorithm does not successfully solve, Value3 is the number of cases where the algorithm could not reach the termination criteria thoughsufficient computational time.TABLE VE
XPERIMENTAL RESULTS IN D ATA S ET UNNING TIME . Size The traditional criteria ( (cid:15) ) The proposed criteria ( K )( − ) ( − ) ( − ) ( − ) ( − ) 5 10 15 20 3010 0.0003 0.0013 0.0045 0.0078 0.0111 0.0002 0.0003 0.0003 0.0004 0.0006100 0.0130 0.0410 0.0913 0.0978 0.1061 0.0043 0.0065 0.0084 0.0105 0.0141500 0.3203 0.9352 1.4508 1.7348 2.4032 0.1049 0.1441 0.1806 0.2179 0.2877 a The values presented in the table are the average values of running time (seconds) of only the success cases from the 50 cases.TABLE VIA
VERAGE RUNNING TIME OF TESTED ALGORITHMS IN D ATA S ET Instacnes Running time (seconds)OPPA EHPA APS OPPA-D (Ours)Da-Com-1 0.0006 0.0025 0.0029 0.0006Da-Com-2 0.0068 0.0038 0.0065 0.0046Da-Com-3 0.0017 0.0040 0.0346 0.0020Da-Com-4 0.1194 0.2675 0.2061 0.1213Da-Com-5 0.2778 0.6585 0.5086 0.2833Da-Com-6 0.4198 1.1887 1.0730 0.4267Da-Com-7 1.5913 11.5127 15.3770 1.5929Da-Com-8 15.2210 216.5146 262.6555 15.3558 a Values presented in the table are the average value of 30 trials. i. When generating the data set, calculate the shortest pathof each graph.ii. Evaluate the algorithm in the data set. In each iteration ofthe algorithm, calculate the L D − P ath . If the L D − P ath ofthis iteration is equal to the shortest path length, recordthe corresponding running time and proceed to Step iii;otherwise, repeat Step ii.iii. Execute the algorithm for 50 or more iterations to ensurethe convergence of the algorithm. If the L D − P ath main-tains the value of the shortest path length for the wholeevaluation period, one could deduce that the T-Point hasalready occurred. Thus, the recorded running time in Stepii is the running time used for comparison; otherwise,return to Step ii.
2) Comparing different accelerated PPAs using the pro-posed method:
Data Set 3:
The instances ‘Da-Com-1’ to ‘Da-Com-8’ men-tioned in Data Set 1 are adopted.
Tested algorithms:
The OPPA, the OPPA-D, the EHPA andthe APS are evaluated. When tested, each algorithm would berun 30 times in each graph.Tables VI and VII presents the experimental results. It is nota surprise that the running time and the number of iterationsof the OPPA-D and the OPPA are approximately the same, asthe difference between the OPPA and the OPPA-D lies onlyin the termination criterion. The similar performance of theOPPA and the OPPA-D under the proposed evaluation methodproves that the method developed by us could exclude the
TABLE VIIA
VERAGE NUMBER OF ITERATIONS OF TESTED ALGORITHMS IN D ATA S ET Instacnes Number of iterationsOPPA EHPA APS OPPA-D (Ours)Da-Com-1 2 2 16 2Da-Com-2 9 9 16 9Da-Com-3 1 1 16 1Da-Com-4 15 15 16 15Da-Com-5 15 15 16 15Da-Com-6 12 12 16 12Da-Com-7 9 9 16 9Da-Com-8 13 13 16 13 a Values presented in the table are the average value of 30 trials. effects brought by different termination criteria.Though the EHPA share the same number of iterations withthe OPPA and the OPPA-D, the running time of EHPA isapproximately 14 times longer than that of the OPPA and theOPPA-D in instance ‘Da-Com-8’. The above findings indicatethat when the effects of the termination criterion is excluded,the EHPA fails to accelerate the PPA as expected; instead, therunning time per iteration is significantly longer.According to Table VI, on average, the running time ofthe APS is even longer than the EHPA. Furthermore, it issurprising that, in all instances, the number of the iterationsof the APS is the same, which is 16. We carefully examinedthe codes provided by the authors [33] and found out that at thebeginning of the iterative procedure of the APS, the algorithmhas a 15-iteration initialization run. One could easily deducethat the APS only uses one iteration to reach the T-Point afterthe initialization. Hence, it suggests that the initialization of theAPS costs the majority of the running time of the algorithm,and further leads to the inefficiency of the APS.In general, the proposed methods do have the ability toavoid the effects brought by different stopping criteria, as theperformance of the OPPA and the OPPA-D is approximatelythe same in the experiment. The most important finding in thisexperiment is that the accelerated OPPAs, i.e., the EHPA andthe APS fail to accelerate the PPA as expected. A potentialreason might be the failure of the authors to not consider thepossible effects brought about by the termination criteria.
EEE TRANSACTIONS 9
V. D
ISCUSSIONS & C
ONCLUSIONS
Though many studies have proposed various solutions toaccelerate the PPA, few have considered the core question:when does the physarum solver distinguish the shortest pathfrom other paths? By defining the concepts of the T-Pointand the D-Path in Section III, we showed that the physarumsolver distinguishes the shortest path at the T-Point. We alsopropose the use of the T-Point as the inherent limit, beyondwhich the algorithm can terminate, hence improving efficiencywhile maintaining accuracy. The mathematical proof remainsin the pipelines for further work. We hope that in the future,the phenomenon of T-Point could also be experimentallyconfirmed in biological experiments.Based on the T-Point and the D-Path, a new terminationcriterion for the OPPA and a corresponding algorithm, OPPA-D, are developed. The OPPA-D is actually the OPPA with thenewly proposed termination criterion. From our experimentsin Section IV-A, the OPPA-D outperforms two state-of-the-art and widely used accelerated OPPAs. In order words, theOPPA significantly outperforms the accelerated OPPAs despitea slight modification to the stopping criterion. This begsanother question: does the accelerated OPPAs beat the OPPAas claimed?In section IV-B, we found that although the proposedtermination criterion is superior to the traditional one, bothof them are sensitive to the predefined parameters, whichwould further affect the validity of comparing different ac-celerated OPPAs. Hence, in order to evaluate the efficiencyof accelerated OPPAs, the effects of the termination criterianeed to be excluded. Thus, we proposed a new evaluationmethod to compare the accelerated OPPAs more objectivelyin Section IV-C1. From our results in section IV-C2 one findsthat the accelerated OPPAs is not as efficient as the OPPA,as previously reported. Previous research did not evaluate thereal performance of the accelerated OPPA by neglecting thepossible effects brought by the termination criteria. However,with the proposed T-Point and the D-Path, and the developedevaluation method, experiments in our studies demonstrate thatthe OPPA is the fastest algorithm among the tested algorithms.We hope the findings in this paper would change the standingcritique that OPPAs are not efficient.A
CKNOWLEDGMENT
The work is partially supported by National Natural ScienceFoundation of China (Grant No. 61973332).R
EFERENCES[1] R. Storn and K. Price, “Differential evolution–a simple and efficientheuristic for global optimization over continuous spaces,”
Journal ofglobal optimization , vol. 11, no. 4, pp. 341–359, 1997.[2] S. Mirjalili, S. M. Mirjalili, and A. Lewis, “Grey wolf optimizer,”
Advances in engineering software , vol. 69, pp. 46–61, 2014.[3] D. Karaboga and B. Basturk, “A powerful and efficient algorithm fornumerical function optimization: artificial bee colony (abc) algorithm,”
Journal of global optimization , vol. 39, no. 3, pp. 459–471, 2007.[4] J. H. Holland, “Genetic algorithms,”
Scientific american , vol. 267, no. 1,pp. 66–73, 1992.[5] M. Dorigo, M. Birattari, and T. Stutzle, “Ant colony optimization,”
IEEEcomputational intelligence magazine , vol. 1, no. 4, pp. 28–39, 2006. [6] S. Mirjalili and A. Lewis, “The whale optimization algorithm,”
Advancesin engineering software , vol. 95, pp. 51–67, 2016.[7] L. Zhang, H. Pan, Y. Su, X. Zhang, and Y. Niu, “A mixed representation-based multiobjective evolutionary algorithm for overlapping communitydetection,”
IEEE Transactions on Cybernetics , vol. 47, no. 9, pp. 2703–2716, 2017.[8] Y. Tian, X. Zheng, X. Zhang, and Y. Jin, “Efficient large-scale multi-objective optimization based on a competitive swarm optimizer,”
IEEETransactions on Cybernetics , vol. 50, no. 8, pp. 3696–3708, 2020.[9] Y. Huang, Y. Gao, Y. Gan, and M. Ye, “A new financial data forecastingmodel using genetic algorithm and long short-term memory network,”
Neurocomputing , 2020.[10] A. W. Mohamed, A. A. Hadi, and K. M. Jambi, “Novel mutation strategyfor enhancing shade and lshade algorithms for global numerical opti-mization,”
Swarm and Evolutionary Computation , vol. 50, p. 100455,2019.[11] P. Hu, J.-S. Pan, and S.-C. Chu, “Improved binary grey wolf optimizerand its application for feature selection,”
Knowledge-Based Systems , p.105746, 2020.[12] M. Wang and H. Chen, “Chaotic multi-swarm whale optimizer boostedsupport vector machine for medical diagnosis,”
Applied Soft Computing ,vol. 88, p. 105946, 2020.[13] A. Tero, S. Takagi, T. Saigusa, K. Ito, D. P. Bebber, M. D. Fricker,K. Yumiki, R. Kobayashi, and T. Nakagaki, “Rules for biologicallyinspired adaptive network design,”
Science , vol. 327, no. 5964, pp. 439–442, 2010.[14] T. Nakagaki, M. Iima, T. Ueda, Y. Nishiura, T. Saigusa, A. Tero,R. Kobayashi, and K. Showalter, “Minimum-risk path finding by anadaptive amoebal network,”
Physical review letters , vol. 99, no. 6, p.068104, 2007.[15] T. Nakagaki, H. Yamada, and ´A. T´oth, “Maze-solving by an amoeboidorganism,”
Nature , vol. 407, no. 6803, pp. 470–470, 2000.[16] T. Nakagaki, “Smart behavior of true slime mold in a labyrinth,”
Research in Microbiology , vol. 152, no. 9, pp. 767–770, 2001.[17] T. Nakagaki, H. Yamada, and A. Toth, “Path finding by tube morpho-genesis in an amoeboid organism,”
Biophysical chemistry , vol. 92, no.1-2, pp. 47–52, 2001.[18] X. Zhang, Q. Wang, A. Adamatzky, F. T. Chan, S. Mahadevan, andY. Deng, “An improved physarum polycephalum algorithm for theshortest path problem,”
The Scientific World Journal , vol. 2014, 2014.[19] V. Bonifaci, “Physarum can compute shortest paths: A short proof,”
Information Processing Letters , vol. 113, no. 1-2, pp. 4–7, 2013.[20] L. Becchetti, V. Bonifaci, M. Dirnberger, A. Karrenbauer, andK. Mehlhorn, “Physarum can compute shortest paths: Convergenceproofs and complexity bounds,” in
International Colloquium on Au-tomata, Languages, and Programming . Springer, 2013, pp. 472–483.[21] A. Karrenbauer, P. Kolev, and K. Mehlhorn, “Convergence of the non-uniform physarum dynamics,”
Theoretical Computer Science , 2020.[22] Y. Liu, C. Gao, Z. Zhang, Y. Lu, S. Chen, M. Liang, and L. Tao, “Solvingnp-hard problems with physarum-based ant colony system,”
IEEE/ACMtransactions on computational biology and bioinformatics , vol. 14, no. 1,pp. 108–120, 2015.[23] X. Chen, Y. Liu, X. Li, Z. Wang, S. Wang, and C. Gao, “A newevolutionary multiobjective model for traveling salesman problem,”
IEEE Access , vol. 7, pp. 66 964–66 979, 2019.[24] Y. Sun, P. N. Hameed, K. Verspoor, and S. Halgamuge, “A physarum-inspired prize-collecting steiner tree approach to identify subnetworksfor drug repositioning,”
BMC systems biology , vol. 10, no. 5, p. 128,2016.[25] S. Xu, W. Jiang, X. Deng, and Y. Shou, “A modified physarum-inspiredmodel for the user equilibrium traffic assignment problem,”
AppliedMathematical Modelling , vol. 55, pp. 340–353, 2018.[26] M.-A. I. Tsompanas, G. C. Sirakoulis, and A. I. Adamatzky, “Evolvingtransport networks with cellular automata models inspired by slimemould,”
IEEE Transactions on Cybernetics , vol. 45, no. 9, pp. 1887–1899, 2015.[27] X. Zhang and S. Mahadevan, “A bio-inspired approach to traffic networkequilibrium assignment problem,”
IEEE Transactions on Cybernetics ,vol. 48, no. 4, pp. 1304–1315, 2018.[28] H. Jiang, X. Liu, S. Xiao, C. Tang, and W. Chen, “Physarum-inspiredautonomous optimized routing protocol for coal mine manet,”
WirelessCommunications and Mobile Computing , vol. 2020, 2020.[29] Y. Song, L. Liu, H. Ma, and A. V. Vasilakos, “Physarum optimization: anew heuristic algorithm to minimal exposure problem,” in
Proceedingsof the 18th annual international conference on Mobile computing andnetworking , 2012, pp. 419–422.
EEE TRANSACTIONS 10 [30] X. Zhang, Y. Zhang, Z. Zhang, S. Mahadevan, A. Adamatzky, andY. Deng, “Rapid physarum algorithm for shortest path problem,”
AppliedSoft Computing , vol. 23, pp. 19–26, 2014.[31] Q. Wang, X. Lu, X. Zhang, Y. Deng, and C. Xiao, “An anticipationmechanism for the shortest path problem based on physarum poly-cephalum,”
International Journal of General Systems , vol. 44, no. 3,pp. 326–340, 2015.[32] Q. Cai and Y. Deng, “A fast bayesian iterative rule in amoeba algorithm.”
International Journal of Unconventional Computing , vol. 14, 2019.[33] C. Gao, X. Zhang, Z. Yue, and D. Wei, “An accelerated physarum solverfor network optimization,”
IEEE Transactions on Cybernetics , vol. 50,no. 2, pp. 765–776, 2018.[34] X. Zhang, A. Adamatzky, X.-S. Yang, H. Yang, S. Mahadevan, andY. Deng, “A physarum-inspired approach to supply chain networkdesign,”