extfieldsofstudy:
Mathematics
Computer Science

text:
Scoring Alternative Forecast Distributions: Completing the Kullback Distance Complex

We develop two surprising new results regarding the use of proper scoring rules for evaluating the predictive quality of two alternative sequential forecast distributions. Both of the proponents prefer to be awarded a score derived from the other's distribution rather than a score awarded on the basis of their own. A Pareto optimal exchange of their scoring outcomes provides the basis for a comparison of forecast quality that is preferred by both forecasters, and also evades a feature of arbitrariness inherent in using the forecasters' own achieved scores. The well-known Kullback divergence, used as a measure of information, is evaluated via the entropies in the two forecast distributions and the two cross-entropies between them. We show that Kullback's symmetric measure needs to be appended by three component measures if it is to characterise completely the information content of the two asserted probability forecasts. Two of these do not involve entropies at all. The resulting 'Kullback complex' supported by the 4-dimensional measure is isomorphic to an equivalent vector measure generated by the forecasters' expectations of their scores, each for one's own score and for the other's score. We foreshadow the results of a sophisticated application of the Pareto relative scoring procedure for actual sequentional observations, and we propose a standard format for evaluation.

Introduction
Bruno de Finetti's "partly-baked" musings with I. J. Good (1962) on the use of a penalty function to gauge the relative predictive qualities of competing forecast distributions have developed by now into a complete theory of proper scoring rules, replete with practical computational procedures. Gneiting and Raftery (2007) present a fairly comprehensive technical review of both the theory and its application. Lad (1996, Chapter 6) presents a pedagogical introduction. Proper scoring procedures are meant to provide the basis for a subjectivist understanding of forecast evaluations that can replace completely the common practice of "hypothesis testing." Such testing is regarded by subjectivists to be without foundation. The very notion that data are generated by an unobservable probability distribution is rejected in their characterisation of probability and its connection with empirical statistical research. In this present tribute to the imagination of de Finetti we propose a striking development in scoring theory. Our analysis gives rise to a completed understanding of Kullback's information measure in the context of the total logl scoring function. Implications for statistical practice are substantial.
The logarithmic scoring function for a probability mass vector (pmv) p N asserted by a forecaster "p" for a quantity X on the basis of an observation X = x o within its realm of possibility R(X ) = {x 1 , x 2 , ..., x N } is specified as the logarithm of the probability assessed for the outcome x o that is actually observed: S(X = x o , p N ) ≡ log(p o ). It has been known since early developments of proper scoring rules by Savage (1971) and by Matheson and Winkler (1976) that a pmv proponent's expectation of one's own log score equals the negentropy value embedded in that pmv. Symbolically, E p [log P(X = x o ))] = ∑ p i log p i . Moreover, p's expectation of the amount by which p's own log score will exceed the score of someone else's assertions, q N , is equal to the relative entropy in p N with respect to q N . That is, This relative entropy measurement is also widely known as the Kullback-Leibler (1951) directed divergence between the pmv's p N and q N , and is commonly denoted by D(p N ||q N ). A non-symmetric function of the two Subsequent consideration of scoring rules by Lad, Sanfilippo and Agrò (2012) has addressed the fact that the log score for mass functions, while proper, is incomplete in that it does not make use of all the observations with their asserted probabilities that are entailed in the assertion of p N and the observation of X . Augmented by the sum of (1 − p i ) log(1 − p i ) values for possibilities of X that are observed not to obtain, the log scoring rule becomes the total log scoring rule. Its prevision (expectation) according to the proponent who asserts p N equals the negentropy plus the negextropy of the asserted mass function. Negextropy was characterized in an article by Lad, Sanfilippo, and Agrò (2015) as a complementary dual of negentropy. The structure of the Kullback-Leibler divergence entertains a complementary dual as well: Parallel considerations define the complementary symmetric Kullback distance as D c (p||q) ≡ D c (p||q) + D c (q||p). Together, the Kullback distance and its complementary distance sum to the Kullback total symmetric distance: The present article presents two new developments in our understanding of associated issues. In the first place, we consider the expectations of the asserters of both p N and q N , each for one's own total log score and for the other's score. Attention turns to these proponents' considerations of the net gains they each might achieve on the basis of their professed willingness to trade their achieved scores for their personally assessed prevision values. (This constitutes the very definition of their previsions.) We find an appropriate strategy for assessing the relative merits of the two pmv assertions: perhaps surprisingly, it is to award to each forecast proponent the net gain achieved by the other's pmv, scaled appropriately, rather than to award to each forecaster one's own achieved proper score.
Secondly, our analysis turns to the pair of expected net gains and expected comparative gains of the forecasters. We learn that Kullback's total symmetric distance measure can be generated equivalently by two other pairs of functions that are different from that which specifies the Kullback-Leibler divergence. This result suggests that Kullback's symmetric distance function should be extended to a four-dimensional vector function if it is to characterise all the distinct features of the information expressed in the two pmv's. It turns out that the vector space spanned by this 4-D vector function is isomorphic to the space spanned by the two proponents' previsions for their own scores and for each other's score.
Section 2 of this article identifies the two pmv proponents' expected net gains and expected comparative gains according to their uncertain assessments of an unknown quantity X . Section 3 develops the logic of the Pareto-optimal assessment of their comparative forecast quality we propose here, while the scale of this comparative scoring procedure is discussed in Section 4. Scaling requires consideration of the variances assessed by the two forecasters for various score components. Attention turns to the higher dimensions of the Kullback complex that arise from these considerations in Section 5, and Section 6 displays why the appended dimensions of the complex are not Bregman divergences. Section 7 portrays the historical context of Kullback's analysis of the problem, and outlines its reconsideration according to the subjectivist understanding. Section 8 presents briefly some experience with and some proposals for application to statistical problems.
Deliberations in this article make use of the syntax preferred for Bruno de Finetti's operational subjective construction of coherent prevision, discussed by de Finetti (1967de Finetti ( , 1972 and Lad (1996, pp. 39-43). Events are characterised as numerical quantities (numbers, not sets) with the realm of possibility R(E) = {0, 1}. Otherwise they are addressed with the same syntax as any other quantity whose realm is finite and discrete. Prevision, P(X ), is a personal value assertion about X , a price at which its proponent is avowedly indifferent to a prize of X or of P(X ). Since a prevision is an assertion made by someone, the P symbol shall be subscripted by a mark denoting who is making the assertion, as in P p ( . ) and P q ( . ) used here. We shall refer to the proponents of these previsions merely as "p" and "q", respectively. One's prevision for an event is also called a probability. Parentheses around an expression that may be true or false denotes an event whose numerical value equals 1 if the expression is true, and 0 if it is false. Thus, for example, E i ≡ (X = x i ) is an event that equals 1 if and only if it happens that X is in fact found to equal x i ; otherwise E i = 0. In the more commonly used measure-theoretic syntax of formalist probability ala Kolmogorov, prevision is equivalent to expectation; in these terms, the expectation of a 0 − 1 random variable is a probability.

Interpersonal assessment of the Total Log Scores
The purpose of this Section is to formalise the notation for the total logarithmic score that will be used to assess our two forecasting distributions (pmvs), and to introduce the formalities of two concepts that will be central to our considerations. The first is the net gain achieved by a forecaster who asserts a pmv for an observable quantity. It defined as the difference between the score that will be achieved and the forecaster's expectation for the score. The second is the comparative gain achieved by one of the forecaster's score relative to that achieved by the other. It is defined simply as the difference between their two scores. We can then identify each forecaster's professed expectation for each of these gains. In doing so we will have identified four fundamental previsions that are relevant to evaluating the forecast quality of the two forecast pmv's on the basis of observations.
Consider the standard problem of statistical forecasting an unknown quantity X whose realm of possible observation values is R(X ) = {x 1 , x 2 , ..., x N }. An associated vector of partition events, E N , is composed of its logically related components E i ≡ (X = x i ) for i = 1, 2, ..., N. The sum of these events must equal 1 since one and only one of them must equal 1. Because the considered value of X is unknown when the forecasters assess probabilities, it is just not known which of the events is the 1 and which of them are the 0's. This will be identified when we observe the value of X . We shall study the uncertain opinions of two people who publicly assert distinct probability mass functions via the vectors p N and q N , respectively. Once the value of X is observed, their assertions will be assessed via a widely touted proper score, the Total Log Scoring function, defined for p N (and similarly for q N ) by as long as the pmv p N is strictly positive and N > 1. Here p o denotes the probability asserted for the value of (X = x o ) that eventually comes to be observed; and the sum in the second line of (1) includes all summands of the The "total log score" is a negative valued score, but it constitutes an award. A larger score (smaller absolute value) is achieved by a more informative distribution. As extremes, the pmv that accords probability 1 to the value of X that occurs, (X = x o ), receives a score of 0; those which accord probability 1 to any particular value of X that does not occur receive −∞. The uniform pmv receives log(1/N) N.B. From this point on we shall no longer write subscripts on vector quantities to denote their size, nor the limits on the range of summations, always understood to be summations over units from 1 to N.
To foreshadow a relevant consideration that will be addressed at the very end of the analysis, we note here that each of the two summands constituting the total log score, ∑ p log p and ∑(1 − p) log(1 − p) is itself a proper score of the pmv p, as is any positive linear combination of these two. Whereas the constructions we now discuss are developed in terms of the total log score, they pertain to each of these addends separately in the same way.

Previsions for one's own score and for another's score
In the following discussion we shall denote the entropy and extropy functions, respectively, by Sums such as − ∑ p log(q) and − ∑ q log(p) are called "cross entropies", and will be denoted by CH(p, q) and CH(q, p). Similarly, "cross extropies" will be denoted by CJ(p, q) ≡ − ∑(1 − p) log(1 − q). Recognise that cross entropies/extropies are not symmetric in their arguments. For example, CH(p, q) = CH(q, p).
Theorem 1. The linearity of coherent prevision identifies the previsions of both p and q for their own scores and for the other's score as: The proof is immediate from the application of linear coherent previsions to the first line of equation (1).
These four previsions exhaust the content of what p and q assert about their own prospective score values and for the scores to be achieved by one another, before the value of X is observed. We shall designate them as the four fundamental previsions, and we shall learn why this nomenclature is appropriate. Let's think a moment about their relative sizes.
Since the scoring function S T L ( . , . ) is proper, each of p and q expects to achieve a greater score than will the other; for the proponent of a pmv expects to achieve a maximum score by publicly professing one's actual pmv, p, rather than any other pmv, q. For this reason, the use of proper scoring rules is said to promote honesty and accuracy in the profession of probabilities. Thus, for example, −[H(p) + J(p)] surely exceeds −[CH(p, q) + CJ(p, q)], a statement which is true when p and q exchange places as well. However, the relative sizes of the two assessors' own expected scores, P p [S T L (X , p)] and P q [S T L (X , q)], may be in any order. This just depends on what they each assert about X , which is p or q, respectively.
When cross-score previsions are considered, p's prevision for q's score may be either smaller or larger than q's expectation for q's own score. In a computational survey of pmv's within the unit-simplex we found that the size of P p [S(X , q)] exceeds P q [S(X , q)] in about 10 to 15 percent of paired (p 3 , q 3 ) selected uniformly at random from the unit simplex. We have not yet characterised the situations in which this occurs, though it should be straightforward.
N.B. From this point we no longer use bold print to distinguish vectors from their components. These will be identifiable by their context.
We are now ready to analyse important features of p's and q's expectations regarding one another's scoring performance.

Expected Net Gains and Expected Comparative Scores
Consider first the "net gains" to be achieved by the proponents of p and q as a result of their receiving the total log scores, S(X , p) or S(X , q) respectively, as a gain in return for paying out the prices at which each of them values this score: (2) Both p and q expect a personal net gain of 0, since their assertions are presumed to be coherent. One's prevised (expected) score is the price at which a person values this unknown score to be achieved. The proponent of a pmv is avowedly willing both to buy and to sell a claim to the score for this price. For example, p is avowedly indifferent to the values of S(X , p) and However, the proponents of p and q do not expect each others' net gains to equal 0. Before assessing them, it is worthwhile to recognize that there is something arbitrary about our definition of Net Gain. For the forecaster is also indifferent to an award of the negative net gain, the result of an exchange in which the forecaster offers to sell the value of S(X , p) in return for receiving P p [S(X , p)]. This would reverse the sign of an awarded Net Gain as we have defined it. This is the sense in which an award to p of a proper score value or its prevision assessment would be arbitrary. Nonetheless, we are now ready to study it as defined.
Theorem 2. P p [NG(X , q)] may be positive or negative valued, as may P q [NG(X , p)].
As proof, note specifically that Now we have already noticed and remarked that p's prevision for q's score may be greater than or less than q's prevision for the same. If p expects a score for q greater than q expects, then p will expect a positive net gain for q because p is thinking that q expects too low a score. In such a case, If p expects a score for q less than q expects, the reverse is be true: The forecasters' expectations of their comparative gains relative to the other's are another matter. Defining the "comparative gain of p over q" on the basis of observing X as it is easy to see that both forecasters expect a positive comparative gain.
Theorem 3. Both forecasters always expect a positive comparative gain for themselves.
As proof, it is clear that and similarly, using the four fundamental previsions we identified in Theorem 1. Both of these previsions are positive because the scoring rule is proper. (3) and (5), it is apparent that summing each of the paired expected net gains and comparative gains yields precisely the negative results of one another. It is worth stating this too as a theorem.

Now comparing the equation pairs numbered
Notice however that the paired summands for these opposite magnitudes are different.
We shall now continue this analysis by devising an improved procedure for the method of comparative scoring, and then we shall use this result to extend our understanding of Kullback's symmetric divergence.

Evaluating pmv's by trading inhering Net Gains
Since p assesses the personal net gain NG(X , p) with prevision 0 while valuing NG(X , q) with the non-zero prevision specified in equation (3), p would surely be willing to trade a claim to plus or minus NG(X , p) in exchange for a claim to NG(X , q) if P p [NG(X , q)] is positive, and to −NG(X , q) if that prevision is negative. For p assesses NG(X , p) with value 0, whereas either NG(X , q) or −NG(X , q) is assessed with positive value. In the same way, q values claims to NG(X , q) with prevision 0, while valuing the net gain of p with the prevision P q [NG(X , p)], as specified in in the lower half of equation (3). Thus both p and q would be pleased to make an exchange in which p receives (appropriately either the positive or negative) value of the NG(X , q) from q while q receives similarly the value of the net gain NG(X , p) (or its negative) from p. In this exchange, both p and q would each be providing the other with something personally regarded as worthless. Thus, both of them would be pleased by the exchanges. Altogether they constitute a "Pareto optimal" exchange.
To clarify this waffle about receiving "a net gain or its negative", note for example that q is avowedly willing to buy or to sell the value of S(X , q) for P q [S(X , q)], since q values these equally. Now with these two options available from q as proclaimed, p would willingly take up q's offer to to sell S(X , q) at this price if P p [S(X , q)] > P q [S(X , q)], or to take up q's offer to buy S(X , q) at this price if P p [S(X , q)] < P q [S(X , q)]. The same structure of relative preference applies to q in reacting to p's proclaimed indifference to S(X , p) and P p [S(X , p)]. To simplify the formalities of this waffling about the action whether to buy or to sell, during further discussion we shall denote the content of such exchanges by saying NG * (X , q) is exchanged for NG * (X , p). This may mean that p gives NG(X , p) to q in return for P p [NG(X , p)] or that p gives P p [NG(X , p)] to q in exchange for NG(X , p), at the choice of q. Simultaneously, q gives to p the value of NG(X , q) in return for P q [NG(X , q)] or vice versa at the choice of p. Of course this specificity would need to be recognized and recorded in the computation of the resulting net gains.
One difference in the assessed values of the two sides of such an exchange is that both p and q would assess their own net gains being given up with a different variance than the net gains they would receive. This constitutes no problem for the usual characterisation of coherent prevision nor for the acceptability of the proposed exchanges. In defining prevision, the scales of exchange are overtly kept small enough that linearity of utility applies to practical valuations. However, we can rescale these agreeable exchanges so that both p and q would regard the two sides of their exchanges symmetrically with the same variance, that is, the same risk. To facilitate this redesign, the pmv proponents would want to rescale the two sides of the exchanges they offer by their assessed standard deviations for them. That is, p would offer to give q NG * (x,p) Both p and q regard both of these offers with positive value, so both would be accepted.
Thus, we have found an intriguing and novel solution to the question of how to evaluate the relative quality of p's and q's pmv assertions when the value of X is actually observed. Simplifying the net result of those two exchanges that are agreeable to both of them, the appropriate way to score the pmv assertions of p and q would be to award to q the score of ,p)] , and to award to p the score of These exchanges constitute for us a Pareto optimal scoring procedure for comparing the relative quality of p's and q's forecasting distributions : to award to q the amount of p's achieved Net Gain * scaled by the sum of the square roots of p's and q's precisions for that Net Gain * , and to award to p the amount of q's Net Gain * scaled by the sum of the square roots of p's and q's precisions for that Net Gain * .
In order to compute these relative information scores, we need to compute four relative variances: Let us digress momentarily to think about these variances in a brief Section of its own.

... and as to the variances
Here are some considerations regarding the cross assessments of the standard deviations. Firstly, remember that NG(X , p)] is a specified number. Now think of the total log score in the form of S(X , p) = log( p o 1−p o ) + ∑ log(1 − p) identified in the third line of equation (1). In this form it is clear that These follow for the same reason, that the summations ∑ log(1 − p) and ∑ log(1 − q) are specified numbers.
according to the standard result for any quantity, that V (Z) = P(Z 2 ) − [P(Z)] 2 .

Kullback's Symmetric Divergence Complex
In this Section we shall learn how and why the structure of Kullback's symmetric divergence function naturally suggests that it be embellished by vectorial components in three independent dimensions.
To begin, recall the details of the Kullback-Leibler directed divergence function D( . || . ) and its extropic complement D c ( . || . ). These functions are also known as the relative entropy and extropy functions of p with respect to q: and similarly for the "reverse directions", Those strict exceedance inequalities presume only that the pmvs p, q are not identical: p = q. Kullback's (1959, p.6) symmetric divergence function D(p, q) is generated by summing the non-symmetric directed divergences: (q||p). Similarly for the complementary symmetric divergence, D c (p, q) ≡ D c (p||q) + D c (q||p). The "total symmetric divergence" is denoted by D T (p, q) = D(p, q) + D c (p, q) .

Expected comparative gains sum to the total divergence
Now an intriguing idea derives from examining once again the expected comparative gain equations, this time algebraically in the context of the Total Log Scoring function, and recognising their relations to the Kullback-Leibler directed divergence functions. We state this as a Corollary 1 to Theorem 3. In the context of the Total Log Scoring function S T L (X , p) introduced in equation (1), the previsions for the comparative gain reduce algebraically by direct substitutions to Each pmv proponent's expected comparative gain over the other's is equal to the relative entropy plus the relative extropy of the proponent's own pmv with respect to the pmv of the other. This is the Kullback-Leibler directed distance, D( . || . ), plus its complementary dual, D c ( . || . ), both of which are known to be positive. So both the proponents of p and q expect a positive comparative gain in total score relative to the other. The positivity of each's expected gain over the other does not depend on which of them expects to receive a larger score personally. Again, both of them expect to receive a larger score than the other (on account of the propriety of the scoring rule) but surely the size of only one of their own actual scores will be larger than that of the other. In fact, their actual comparative gains are the negative equivalents of one another, according to their defining equation (4).
It is now worth noting that the sum of p's and q's previsions for their comparative gains can be seen in equations (5) to equal the total symmetric divergence of the pmv's.
Corollary 2 to Theorem 3. Again in the case of the Total Log Score, This is despite the fact that the sum of their actual comparative gains equals 0. Further investigation yields still another intriguing equivalence.

Expected Net Gains sum to the same thing
Intrigue arises from examining now the algebraic details of p's and q's expectations of each others' net gains, viz., , and similarly Each of these is the difference between one's own expection and the other's expectation of the other's log odds for the occuring event.
Now summing these expected net gains yields a companion to Corollary 2 of Theorem 3, evident from basic definitions: Corollary 3 of Theorem 3. Again using the Total Log Scoring rule,

An alternative generator of D T (p||q).
Perhaps surprisingly, at this late stage in the study of the divergence function, we can now recognise that still two more companion functions can be identified whose sums yield the same symmetric divergence as do the Kullback-Leibler functions D(p||q) and D(q||p). The awareness arises from Corollaries 2 and 3 of Theorem 3. The first pair, also directed functions, arises formally from switching the components of the finite discrete entropy function, −H(p) = ∑ p log p that get "differenced" in the production of D( . || . ). Rather than thinking about ∑ p (log p − log q), we consider differencing the "p" instead of the "log p", raising consideration of ∑(p − q) log p. We shall denote this alternative directed function by ∆( . || . ), motivating it below: and then similarly for the "reverse directions", Referring to our formulation of cross prevision assertions, it is apparent that ∆(p||q) = CH(q, p) − H(p), the amount by which q's cross entropy for p exceeds the entropy in p itself. Similarly, ∆ c (q||p) = CJ(q, p) − J(p), the amount by which q's cross extropy for q exceeds q's own extropy.
These considerations yield a surprising result: both the directed distance D(p||q) and its reverse D(q||p), and this new function ∆(p||q) and its reverse ∆(q||p), sum to the same thing, ... Kullback's symmetric divergence, which we are denoting by D(p||q). The same is true of the complementary sums, which we denote by D c (p||q). This result arises as Corollary 1 to Theorem 4.
Here the sums of the ∆( . || . ) functions apply the total log score to the Net Gains, while the D( . || . ) functions arise from its application to the Comparative Gains.

Identifying still a third function that yields D T (p, q)
Expanding the identical sums D(p||q) + D(q||p) and ∆(p||q) + ∆(q||p) which constitute D(p, q) into their component pieces, and then exhibiting them across a line makes it apparent that each of the functions, D(p||q) and ∆(p||q) arises by selecting two of the four summands to define the function, and thus the remaining two to define their "reverses", D(q||p) and ∆(q||p).

D(p||q)
is the sum of the first two pieces of the first line, leaving D(q||p) composed of the second two; while ∆(p||q) is the sum of the first and third piece, leaving ∆(q||p) to be composed of the second and fourth piece. Now seeing these sums and these function pairs in this way suggests that a third (and the accompanying sixth) choice of two of the four summands may be worthy of study as well. Consider using the interior two summands (the second and third pieces of the first line) to constitute still another form of function that would sum (along with the remaining pair) to yield the symmetric divergence D(p||q). This choice would be the sum of the "cross entropy" functions, CH(p||q) +CH(q||p). One difference of this function from D(p||q) and ∆(p||q) is that this sum of cross entropies is already symmetric on its own, whereas the latter two functions require summation with their companions D(q||p) and ∆(q||p) to construct a symmetric function. A second difference of the sum CH(p||q) +CH(q||p) from the D(p||q) and ∆(p||q) functions appears in its location for self-divergence. Whereas the values D(p||p) = ∆(p||p) = 0, the value of CH(p||p) +CH(q||q) = H(p) + H(q), the summed entropies in p and in q. However, defining the symmetric divergence function more generally than heretofore so to eliminate this location shift of "self-distance" yields a pleasing symmetrical result.
Theorem 5. Based on a generalised definition of symmetric divergence, so to account for the location of "self-divergence", Of course those self-divergence values in the first two lines equal 0, whereas in the third line CH(p||p) = H(p) and CH(q||q) = H(q).
The same procedures of deconstruction and reconstruction can be performed with the complementary functions (D c ( . || . ) and ∆ c ( . || . )) in a second line: D c (p||q) is the sum of the first two pieces of this second line, while ∆ c (p||q) is the sum of the first and third pieces.
We propose that the Kullback directed distance function needs to be appended by all three of its three generators if it is to characterise the information content of the statistical problem the forecasters address. Together they constitute a complex.
One might be puzzled. "Who cares about all these generators of the same symmetric distance? The symmetric distance is what it is." The answer derives from an identifiable linear relation between the four functions D T (p, q), D(p||q) + D c (p||q), ∆(p||q) + ∆ c (p||q), CH(p||q) +CJ c (p||q) and the "four fundamental previsions" we introduced in equation (2). We address this now.

The isomorphism of the Kullback complex of generators with the four fundamental previsions
Each of these eight component functions in a Kullback complex is a different linear combination of the various entropies, extropies, cross entropies and cross extropies that have been central to our scoring analysis. These combinations happen to be ordered in such a way that the two spaces spanned by the 4-dimensional vector functions are isomorphic, related by the linear equations What this relation tells us is that each of the three distance function generating pairs measures a distinct information source for understanding of the full information content of Kullback's symmetric distance. The symmetric distance measure is incomplete on its own: it needs to be supplemented by three more components if it is to represent the information content of the pmv assertions p and q. In all, it is identifiable better as a vector which we might refer to as the Kullback complex ... [D(p, q), D(p||q), ∆(p||q),CH(p||q) +CH(q||p)], rather than merely by the symmetric function value D(p, q) alone.
One further recognition is in order. We have couched this discussion in the context of the total logarithmic score (1) and the Kullback total symmetric distance function D T (p||q). It is worth remarking that the total log score can be considered to be the sum of two separate components: the log score itself, ∑ E i log p i , and the complementary log score, ∑Ẽi log (1 − p i ). Actually, both of these components of the total log score are proper scores on their own, as is any linear combination of them. The two component scores assess different aspects of the asserted pmv's p and q. One piece focuses more directly on the center of the mass function, while the other focuses more on the tails. Because these aspects of a pmv may concern different people in different ways with respect to utility, it is usually best to keep track of the two components of the score separately. If the scoring function S T L (X , p N ) were considered to be the sum of the two separate scores, say, S log (X , p N ) + S complog (X , p N ), the isomorphism equations (9) hold separately for each component of the total log score and their previsions too. Tracking these components separately, the total symmetric Kullback complex could be considered to be an 8-dimensional space of measurement functions.
We shall now address a final technical note concerning the relation of these newly recognised dimensions of the information complex to the theory of Bregman divergences.


extfieldsofstudy:
Computer Science
Mathematics

text:
Mapping the Global Structure of TSP Fitness Landscapes

The global structure of combinatorial landscapes is not fully understood, yet it is known to impact the performance of heuristic search methods. We use a socalled local optima network model to characterise and visualise the global structure of travelling salesperson fitness landscapes of different classes, including random and structured real-world instances of realistic size. Our study brings rigour to the characterisation of so-called funnels, and proposes an intensive and effective sampling procedure for extracting the networks. We propose enhanced visualisation techniques, including 3D plots and the incorporation of colour, sizes and widths, to reflect relevant aspects of the search process. This brings an almost tangible new perspective to the landscape and funnel metaphors. Our results reveal a much richer global structure than the suggestion of a ‘big-valley’ structure. Most landscapes of the tested instances have multiple valleys or funnels; and the number, disposition and interaction of these funnels seem to relate to search difficulty on the studied landscapes. We also find that the structured TSP instances feature high levels of neutrality, an observation not previously reported in the literature. We then propose ways of analysing and visualising these neutral landscapes.

by considering the distribution of its local optima. Boese et al (1994) conjectured that the search space of travelling salesman instances under 2-exchange moves has a 'globally convex' or 'big-valley' structure, in which local optima are clustered around one central global optimum. This globally convex structure has subsequently been observed in other combinatorial fitness landscapes such as the NK model (Kauffman and Levin 1987), graph bipartitioning (Merz and Freisleben 1998), and flowshop scheduling (Reeves 1999). Under this view, there are many local optima but they are easy to escape from, with the coarse level gradient leading to the global optimum. This hypothesis has become generally accepted and has inspired the design of local search heuristics referring to a similar principle with different names: adaptive multi-starts, large-step Markov chain, soft-restarts, chained local search and iterated local search. The notion of a big-valley is related to the notion of a 'funnel' structure from the study of energy landscapes in theoretical chemistry, as discussed further in Section 2.2.
However, recent studies on TSP landscapes have revealed a more complex picture (Hains et al 2011;Ochoa and Veerapen 2016b). The big-valley seems to decompose into several sub-valleys or multiple funnels. This helps to explain why certain iterated local search heuristics can quickly find high-quality solutions, but fail to consistently find the global optimum in cases where the global optimum is known. A similar multi-funnel structure has been observed on some continuous optimisation problems (Locatelli 2005;Lunacek and Whitley 2006;Lunacek et al 2008), where its impact on search difficulty has been established. In particular, landscapes with more than one funnel, where the global optimum is located in a deep, narrow funnel are significantly harder. The literature on characterising the multi-funnel structure of combinatorial landscapes is mostly lacking. This is partly due to the lack of adequate tools to study their complexity. We propose using local optima networks to analyse and visualise the global structure of combinatorial fitness landscapes. Local optima networks compress the whole search space into a graph, where nodes are local optima and edges are transitions among them with a given search operator (Tomassini et al 2008;Verel et al 2011). The model emphasises the number, distribution and most importantly, the connectivity pattern of local optima. Modelling landscapes as networks introduces a new set of metrics to analyse fitness landscapes and, interestingly, the possibility of visualising them.
This article complements and extends our recent work on local optima network analysis of TSP fitness landscapes (Ochoa and Veerapen 2016a,b). Our previous work considered only a handful of instances of size up to 700 cities, and deliberately excluded degenerate instances (i.e. instances featuring neutrality). The sampling process was less intensive, and the characterisation of funnels incomplete. The main contributions of this article are: 1. An intensive sampling procedure, which allows the extraction of local optima networks for larger TSP instances (of up to 1500 or so cities) of different types, including both random and structured real-world instances (city instances and drilling problems).
The remainder of this article is organised as follows. The next section gives relevant background on TSP solvers and the notion of 'funnel'. Section 3 defines the local optima network model considered. Section 4 describes the proposed empirical methodology, including the procedures and conceptual tools for extracting the network data, detecting the funnel floors and calculating the funnel basins. It also describes the TSP instances studied. Section 5 overviews our results including an analysis and visualisation of local optima networks, a more traditional fitness-distance analysis, a study of the impact of the sampling parameters and a correlation study between heuristic search performance and global landscape metrics. Finally, Section 6 summarises our main findings and suggests directions of future work.

TSP solvers
This section describes the TSP solvers used in our study both for comparison purposes and as part of the sampling procedure implemented.

Concorde
Concorde is currently the best-performing exact TSP solver (Applegate et al 2007(Applegate et al , 2003a. It has been used to solve the largest non-trivial TSP instances (of up to 85 900 cities) for which provably optimal solutions are known. Concorde is based on a complex Branch & Cut algorithm that uses a multitude of heuristic mechanisms to achieve good performance on a wide range of TSP instances. For example, it carries out a limited number of iterations of the Chained Lin-Kernighan heuristic (described below) during the initial stages of its computation to determine an initial upper bound to the objective value. Additionally, an exact mixed integer program solver is used to compute and refine the lower bound by solving a relaxed linear program of the problem.

The Chained Lin-Kernighan heuristic
The Lin-Kernighan (LK) algorithm (Lin and Kernighan 1973) is a powerful and wellknown heuristic for finding approximate solutions to the TSP. For about two decades, it was the best local search method, and nowadays it is a key component of many state-of-the-art TSP solvers. LK-search is based on the idea of k-exchanges: take the current tour and remove k different links from it, which are then reconnected in a new way to achieve a legal tour. Figure 1a illustrates a 2-exchange move. A tour is considered to be 'k-opt' if no k-exchange exists which decreases its length. LK-search applies 2, 3 and higher-order k-exchanges. The order of a change is not predetermined, rather k is increased until a stopping criterion is met. Thus many kinds of k-exchanges and all 3-exchanges are included. There are many ways to choose the stopping condition and the best implementations are rather involved. We use the implementation available in the Concorde software package (Applegate et al 2003a), which uses do not look bits and candidate lists.
The overall tour-finding strategy using LK-search was, previously, to repeatedly start the basic LK routine from different starting points keeping the best solution found. This practice ended in the 1990s with the seminal work of Martin et al (1992), who proposed the alternative of kicking (perturbing slightly) the LK tour and reapplying the algorithm. If a better tour is produced, the old LK tour is discarded and the new one kept. Otherwise, the search continues with the old tour and kicks it again. This simple yet powerful strategy is nowadays known as iterated local search. It was named Chained Lin-Kernighan (Chained-LK) by Applegate et al (2003b), who also provided an improved implementation to solve large TSP instances. The kick or escape operator in Chained-LK is a type of 4-exchange (depicted in Fig. 1b), named double-bridge by Martin et al (1992). It consists of two 2-exchanges, each of which is a 'bridge' as it takes a legal, connected tour into two disconnected parts. The combination of both bridges must then be chosen in order to produce a legal final tour.

The notion of funnel
The intuition behind the concept of a 'funnel' is captured by Figure 2 where two funnels are depicted as two groups of local optima which are close in configuration space within a group, but well-separated between groups. The term 'funnel' was introduced in the protein folding community to describe "a region of configuration space that can be described in terms of a set of downhill pathways that converge on a single low-energy structure or a set of closely-related low-energy structures" (Doye et al 1999). It has been suggested that the energy landscape of proteins is characterised by a single deep funnel, a feature that underpins their ability to fold to their native state. In contrast, some shorter polymer chains (polypeptides) that misfold are expected to have other funnels that can act as traps. Approaches to elucidate the global landscape structure have led to the concept of disconnectivity graphs (Becker and Karplus 1997;Doye et al 1999;Wales 2005), also known as barrier trees (Flamm et al 2002). In the context of energy landscapes, Berry and Kunz (1995) first introduced the term monotonic sequence to describe a sequence of local minima where the energy of minima is always decreasing. The set of monotonic sequences that lead to a particular minimum was termed 'basin'; in this sense a 'basin' is analogous to a protein folding 'funnel'. The collection of local optima associated to a funnel has also been termed 'super-basin' in the literature introducing disconnectivity graphs (Becker and Karplus 1997).
Energy landscapes in theoretical chemistry and fitness landscapes in optimisation are conceptually related, as has been already observed by Stadler (2002). This relationship is particularly close for continuous optimisation. Locatelli (2005) studied the sources of difficulty in continuous optimisation and finds that it is not strictly related to the number of local optima, but to how chaotic their positions are. Lunacek and Whitley (2006) propose a metric, dispersion, that quantifies the proximity of the best regions in the search space. A high dispersion metric indicates the presence of multiple funnels. In a follow up work, Lunacek et al (2008) studied abstract landscapes with two funnels and find that evolutionary algorithms mostly fail when the global optimum is located in a proportionally smaller funnel. Recent work on exploratory landscape analysis of continuous search spaces reveals that multimodality and global structure are among the most important high-level properties that help differentiate between problem classes (Bischl et al 2012;Kerschke et al 2015).
The literature is much more scarce for discrete search spaces. The notion of a big-valley discussed in the introduction is clearly related to the notion of a single funnel structure, and fitness distance scatter plots and correlation metrics are a standard tool in landscape analysis. Other approaches to understanding the global structure of combinatorial landscapes have been applied to small and simplified problems. Barrier trees have been applied to discrete optimisation problems where the notions of local optima, basins and saddle points are clearly defined, for instance in the context of spin-glasses (Hordijk et al 2003). Flamm et al (2002) have extended these definitions so that barrier trees can be constructed for highly degenerate problems (i.e landscapes with neutrality). They present empirical results for binary strings of up to length 10. Hallam and Prugel-Bennett (2005) construct barrier trees for MAX-SAT problems with up to 40 variables using branch-and-bound to find only the best local optima in the space. Daolio et al (2011) studied the community structure of local optima networks on two classes of instances of the quadratic assignment problem. The two problem classes give rise to different configuration spaces, with the so-called real-like instances revealing a modular structure. The approach is based on a full enumeration of local optima. Therefore, instances of size up to 10 were analysed. In a follow up work with a data-driven approach, the modularity of instances up to size 32 was studied (Iclanzan et al 2014). This work, however, did not relate the community structure to the notion of funnels. Herrmann et al (2016) recently established a connection between groupings (communities) in local optima networks and the notion of funnels. They extracted the networks of NK landscapes of length 20 and various levels of epistasis, and applied a community detection algorithm. Results confirm that landscapes consist of several clusters and the number of clusters increases with the epistasis level. A higher number of clusters, and a larger size of the cluster containing the global optimum were found to lead to a higher search difficulty.

The Local Optima Network Model
This section describes the local optima network model used in our study. We start by defining the notion of fitness landscapes, and follow by formalising the notions of nodes and edges of the network model.
A fitness landscape (Stadler 2002) is a triplet (S, N, f ) where S is a set of potential solutions i.e. a search space; N : S −→ 2 S , a neighbourhood structure, is a function that assigns to every s ∈ S a set of neighbours N(s), and f : S −→ R is an objective function (also called fitness function) that can be pictured as the height of the corresponding solutions.
The search space of a TSP instance of size m is the set of permutations of the m cities. The objective function f is given by the length of the tour, which is to be minimised. In order to model TSP fitness landscapes, we adapted the local optima network model with escape edges (Verel et al 2012). To construct these networks, we need to define their nodes and edges. The definitions are related to the search operators being modelled, specifically, the local search heuristic and escape operators. Our study considers those within the Chained Lin-Kernighan algorithm, namely, the Lin-Kernighan local search heuristic and the double-bridge escape move (described in Section 2.1.2).
Local optima. A local optimum, which in the TSP is a minimum, is a solution s * such that ∀s ∈ N(s * ), f (s * ) ≤ f (s). Notice that the inequality is not strict, in order to allow the treatment of the neutral landscape case.
The neighbourhood N is imposed by LK-search, which considers variable values of k. The local optimality criterion is, therefore, rather stringent. Only a small proportion of all possible solutions are LK-optimal. The set of local optima, which corresponds to the set of nodes in the network model, is denoted by L and its cardinality by n.
Escape edges. Edges are directed and based on the double-bridge operator. There is an escape edge from local optimum x to local optimum y, if y can be obtained after applying a double-bridge kick to y followed by LK-Search. Edges are weighted with estimated transition probabilities between the connected nodes. These probabilities are estimated by the sampling process. Specifically, edge weights are integers indicating the number of times an edge was visited during the sampling process (described in Section 4.1). The set of escape edges is denoted by E.
Local Optima Network (LON). A local optima network is a graph LON = (L, E) where nodes are the local optima L, and edges E are the escape edges. Edges are directed and weighted. Weights indicate transition probabilities.

Empirical Methodology
Our approach extracts and analyses local optima networks of TSP instances of realistic sizes and different types. The objective is to map and characterise the landscapes' global structure. Clearly, a full enumeration of the local optima for TSP instances of non-trivial size becomes unmanageable. Therefore, networks are constructed from a sample of high-quality local optima in the search space. This section starts by describing our sampling methodology and procedure for constructing the LONs. Thereafter, we describe the approach for characterising the funnel structures, which requires both detecting the funnel floors, and computing their basins.

Sampling the network data
To extract the network data, we considered the Chained-LK implementation of Concorde (see Algorithm 1). Instead of discarding the search history, we store, in L, every new improving LK local optima found during the search process. We also create and store, in E, a directed edge between the starting and ending optima after a doublebridge followed by LK. Edges only go from a node with higher cost to a node of equal or lower cost, in order to store the monotonic sequences and thus identify funnels (see Section 4.2). We keep count of the number of times an edge is visited, this number is stored as the edge's weight.
A thousand independent runs of Chained-LK are executed for each TSP instance. We consider two initialisation mechanisms, one producing better initial solutions than the other, in order to have a broader picture of the search space. Half of the runs use the Quick-Borůvka method, the default initialisation for Concorde's Chained-LK, Data: I, TSP instance Result: L, set of local optima, E, set of escape edges Algorithm 1: Local optima network sampling combining 1000 runs of Chained-LK.
which is based on the minimum-weight spanning tree algorithm of Borůvka (Applegate et al 2003b). The other half starts from a random solution. The default kicking procedure in Concorde's Chained-LK is used: the edges involved in the double-bridge are selected using random walks along connected vertices. Each Chained-LK run continues until at least 10 000 consecutive iterations are performed without finding an improving solution. The network is thereafter created by the combination of the unique nodes and edges produced by this sampling process.

Identifying funnel structures
The challenge is to devise an approach for identifying funnel structures, for different types of TSP instances, once the local optima networks have been extracted. Our approach adapts the notion of monotonic sequences from theoretical chemistry (Berry and Kunz 1995). We consider a monotonic sequence as a sequence of local optima where the evaluation of solutions is non-deteriorating. The collection of monotonic sequences leading to the same lowest minimum correspond to the so-called 'monotonic sequence basins' (Wales 2005). These structures have also been called 'superbasins' in the theoretical chemistry literature (Becker and Karplus 1997). We choose here to call them 'funnel basins' or simply 'funnels' borrowing from the protein folding literature. We can distinguish the primary funnel, as the one involving monotonic sequences that terminate at the global optimum. The primary funnel is separated from other neighbouring secondary funnels by transition states laying on a so-called 'primary divide' (Berry and Kunz 1995). Above such a divide, it is possible for a local optima to belong to more than one funnel thorough different monotonic sequences.
Our approach requires us to empirically locate the lowest cost minima which potentially lie at the bottom funnels, and thereafter compute each funnel basin. These procedures are described below.

Identifying the funnel floors
Hains et al (2011) considered funnel floors (or funnel bottoms) as those solutions empirically found after considerable search effort. Specifically, for each TSP instance, Chained-LK was run until at least 10 000 iterations without finding an improving tour, and this procedure was repeated 1000 times from different starting solutions. However, runs were considered separately and the intermediate local optima and transitions among them discarded, therefore, the number of funnel floors was overestimated.
We propose a refinement of this approach to more accurately estimate the number of funnels. We keep a similar computational effort, but store all the different local optima visited across the 1000 runs, and combine them into a single local optima network. This allows us to merge local optima transitions found across different runs. We name the solutions found at the end of each run (i.e. 10 000 consecutive iterations without an improvement) as attractors rather than as funnel bottoms. Our local optima network data shows two interesting features. First, solutions attracting the search process are not always single solutions, but are often part of connected local optima plateaus. Second, attractor nodes are not always at the bottom of funnels: a single run is trapped in an attractor, but when combining the 1000 runs into a single LON, many attractors show escape paths to solutions with lower evaluation. Therefore, in order to refine the process of detecting the conjectured funnel floors, we propose constructing a sub-network from the LON, where nodes are the attractors: the Attractor Network (AN). Once this network is constructed for a given instance, we detect its sinks, i.e. nodes without outgoing edges. We consider these sinks as the solutions at the bottom of funnels, which produced a reduced number of funnels for the studied instances as compared to previous work (Hains et al 2011;Ochoa and Veerapen 2016b). We overview below the notation and procedures used to detect the number of funnel floors.
-Attractor nodes are empirically determined as those local optima at which Chained-LK stalls after a large search effort (10 000 consecutive iterations without finding an improving solution in our implementation). The set of attractor nodes is denoted by A, and its cardinality by a. -A sink node, is an attractor node without outgoing edges to other attractor nodes of lower evaluation. Sink nodes are conjectured to be at the bottom of a funnel structure. The set of sink nodes is denoted by S and its cardinality by s. -A plateau node, is a higher-level node compressing a group of local optima with the same objective function value belonging to a connected component according to the escape edges defined in Section 3. Plateau nodes can also be characterised as attractors or sinks. -Attractor-plateau nodes are calculated by compressing connected attractor nodes at the same objective function level. The set of attractor-plateau nodes is denoted by A p , and its cardinality by a p .
-A sink-plateau node, is an attractor-plateau node without outgoing edges to other attractor-plateau nodes. The set of sink-plateau nodes is denoted by S p and its cardinality by s p .
Our previous local optima network models of TSP landscapes (Ochoa and Veerapen 2016b) involved a less extensive sampling and did not consider plateau nodes. However, the study of a more varied set of instances revealed neutrality (i.e. different solutions with the same objective function value) on the landscapes of most structured instances, especially those of drilling problems. The random instances do not show neutrality, but many of the city instances and a large proportion of the drilling problems do reveal plateaus that can be extensive. It was, therefore, necessary to include plateau nodes, and to define two auxiliary network models (described below) in order to algorithmically detect the sink nodes, and thus funnel floors.
-Attractor Network (AN), the sub-graph AN = (A, E a ) of LON where nodes are the attractors A, and edges E a ⊆ E are the escape edges connecting them. -Attractor-plateau Network (AN p ), the graph AN p = (A p , E p ) formed by contracting the nodes and edges of the Attractors network AN. Figure 3 illustrates the attractor and attractor-plateau networks for the well-studied TSPLIB instance att532 (more details in Table 1), which consists of the 532 largest cities of the USA and considers pseudo-Euclidean distances between them. The procedure by Hains et al (2011), estimated four funnels for this instance. Our analysis, instead, suggests that it has only two funnels, whose sinks are represented in red and yellow in Figure 3. As the figure indicates, there are 48 attractor nodes (a), which are compressed into 7 plateau-attractor nodes (b). Not all the 7 plateau-attractor nodes are funnel floors, as most of them have outgoing edges to other nodes. Instead only two of them, the sinks coloured in yellow and red in plot (b), are conjectured to be at the bottom of funnels. The whole sampled local optima network for this instance, including nodes within 1% in evaluation from the global optimum cost, is visualised in Fig. 5 (a).

Identifying the funnel basins
Once the funnel sinks are detected, we can proceed to identify the funnel basins (see Algorithm 2). This is done by finding all the local optima in the network which are reachable from each funnel sink (sink-plateau for the instances with neutrality). Breadth-First-Search is used for this purpose. The set of unique solutions in the combined paths to a given funnel sink corresponds to the funnel basin. The cardinality of this set corresponds to the funnel size. Notice that the membership of a solution to a funnel might be overlapping, that is, a solution may belong to more than one funnel, in that there are paths from that solution to more than one funnel sink. The relative size of the primary funnel (or any other secondary funnel) is calculated as its size divided by the total number of local optima.
As described above, funnel sinks are nodes without outgoing edges in the attractor network, which is a sub-graph of the sampled local optima network (LON). When considering the whole LON, out-going edges might exist between sub-optimal sinks  (a), and the same network (b) after compressing plateau nodes for a selected instance with 532 cities (att532). In the attractor network, each node is a local optimum with size proportional to its in-strength (i.e. the weighted in-degree). In the attractor-plateau network, nodes are rectangular with width proportional to the number of nodes in the plateau. The instance contains two funnel sink-plateaus, visualised in red (the global optima) and yellow (the secondary funnel sink), the legend indicates their cost values. Grey nodes are above (i.e higher cost) than their respective sinks. and local optima with lower evaluation, which are not sinks themselves, but belong to a different basin. In consequence, before identifying the funnel basins for each sub-optimal sink, we remove all its outgoing edges. This is done by the function disconnectSinks(LON, S) in Algorithm 2. Table 1 summarises the TSP instances studied. We consider 20 instances of moderate size (in the range of 500 to 1500 or so cities) and different types. The first 10 instances are randomly generated using the DIMACS generator 1 . Half of these are composed of uniformly distributed cities (prefixed by 'E'), while in the other half, the cities are clustered (indicated by a 'C'). The suffix number '.x' in the instance name indicates, as per DIMACS convention, a seed of x + 10 000. These synthetic instances are part of a larger set of 200 instances that we generated, with the number of cities being uniformly selected in the range [400, 1600]. These are used in Section 5.6 for a correlation study between landscape metric and heuristic search performance.

Instances
Data: LON: sampled local optima network, S: funnel sinks Result: bsizes: basin sizes vector, basins: funnel basins vector, boverlap: set of local optima in more than one basin Algorithm 2: Identifying funnel basins. Function disconnectSinks(LON, S) removes all outgoing edges from sub-optimal sinks in LON.
The bottom 10 instances in Table 1 are well-known instances from TSPLIB (Reinelt 1991). A popular way of constructing TSP instances is to choose a set of actual cities and define the cost of travel between any two cities as the distance between them. The first 5 TSPLIB instances are constructed in such a way. The last 5 arise from the task of drilling holes in printed circuit boards. The types of edge weights are as follows. EUC-2D refers to the Euclidean distance of points in a 2D plane rounded to the nearest integer. ATT refers to a pseudo-Euclidean distance where the sum of the squares is divided by 10 and the square root of this value is then rounded to an integer. GEO refers to the integer geographical distance computed from latitude and longitude coordinates on the surface of a sphere representing an idealised Earth.
The third column in Table 1 reports the success rate of the 1000 Chained-LK runs used for extracting the network data. By success rate, we mean the ratio of runs that found at least one global optimum. The last two columns give information on the solving difficulty of each instance. Specifically, we report the mean run time and the mean number of branch-and-bound nodes required by Concorde (interfaced with IBM ILOG CPLEX 12.6 as its mixed integer solver) to solve the instances to optimality on a 3.4 GHz Intel Core i7-3770 CPU across 10 runs. Although Concorde is an exact solver, the means are computed since it uses Chained-LK to generate initial solutions. Therefore, times include the Chained-LK time to compute initial solutions. This leads to different execution times and branch-and-bound trees. A number of branch-andbound (B&B) nodes equal to 1 indicates that the lower and upper bounds found in the initial stages of the Concorde solving process match, and thus no actual tree was explored.

Results
The results and analysis are reported in the following six subsections. We start by reporting general statistics of the sampled local optima networks across the instance set. Subsection 5.2 concentrates on features characterising the landscapes' global struc- ture. Subsection 5.3 visualises the local optima network of some selected instances. Thereafter, Subsection 5.4 provides a more classic fitness landscapes analysis. We finish with a study of the effects of sampling parameters on the landscape metrics (Subsection 5.5), and a correlation study between landscape metrics and heuristic search performance (Subsection 5.6). Table 2 reports basic statistics of the sampled local optima networks, including the number of unique global optima go, the number of different local optima n, the number of unique objective function evaluations evals, and the relationship between these two values n/evals as an indication of the amount of neutrality in the landscape. The table also reports the number of connected components c, the proportion of nodes in the network from where there is a path to a global optimum p go ; and the averagel go and maximum l go path length from any node in the network to a global optimum. The last two columns report the averaged and maximum d max in-strength (i.e. incoming weighted degree, where the weight of an edge is the number of times it was traversed during sampling) of nodes in the networks. Results show that, as expected, the number of global optima and the average and maximum path length to a global optimum generally increase with the number of cities for each instance class. There are, however, striking differences among the classes. The randomly generated instances always show a single global optimum and each local optima has a different objective function value (as indicated by columns evals and n/evals). This is not the case for the structured instances, where several global optima are generally the norm. Indeed global optima seem to be located in large plateaus (of several thousands of nodes) for some of the drilling instances. The high amount of neutrality is also revealed by the smaller number of objective function levels as compared to the number of optima, with the ratio being as large as several hundreds or thousands. This is probably because several pairwise distances between cities are the same, and some of these distances are rather small. Therefore, several city orderings will have the same evaluation. The structured instances also show longer path lengths to a global optimum, which can be in part explained by high amount of neutrality.

General network metrics
The average in-strength,d, is usually quite low since most nodes are visited only once. Higher values indicate either that multiple runs find the same nodes or that there is cycling between nodes of equivalent evaluation. This cycling is also why d max is quite high for some instances.
A surprising consequence of sampling and modelling the local optima networks for the very neutral instances such as u1060 and fl1577, is that the number of connected components (c in Table 2) equals the number of runs in the sampling process (1000). This means that each run traverses a different set of solutions, in other words, there are no shared local optima among the runs, even though many runs reach the same objective function values. We suspect that these instances contain such large plateaus, that runs (starting from different initial points) explore completely different parts of them. For example, the number of different global optima found by our sampling process on instance u1060 is 163 569. The metric n/evals in Table 2, gives an estimate of the average size of plateaus for each instance, which can be as large as several hundred. Table 3 reports metrics describing the global structure of the random instances. We report the proportions of: solutions in the primary funnel (i.e the funnel containing the global optimum) f go , solutions in the largest funnel (which can be the primary funnel) f lg , and solutions belonging to more than one funnel f ov . The table also reports the number of attractors a and funnel sinks s. The proportion of the in-strength (i.e. the weighted in-degree) of global optimum sinks to the in-strength of all sinks, d go , is given and the Chained-LK success rate from Table 1 is reproduced in the last column for comparison purposes.

Funnel metrics
Results suggest that the number of funnels generally increases, and the size of the primary funnel decreases, with increases to the number of cities on both instance Table 2: Local optima network metrics. Number of unique global optima go, number of unique local optima n, number of unique objective function evaluations evals, relationship between number of nodes and unique evaluations n/evals, number of connected components c, proportion of nodes from where there is a path to a global optimum p go , average path length to a global optimuml go , maximum path length to a global optimum l go , average in-strength of nodesd and maximum in-strength d max . classes. The random uniform instances (prefixed with 'E'), show a larger number of funnels and a smaller primary funnel, when compared to the clustered instances (prefixed with 'C'). These metrics change faster with the instance size on the uniform instances, specially the number of funnels. The proportion of solutions in more than one funnel is also larger for the uniform instances, which is expected given the larger number of funnels. By virtue of its definition, the relative in-strength of globally optimal sinks is closely related to success rate, but it is not exactly the same, as can be observed from the table. Table 4 reports global metrics for the structured instances. Since these instances contain neutrality, the attractor-plateaus networks (AN p ) are constructed in order to identify the plateau-sink nodes, and thus the number of funnels. Results suggest that the number of funnels is less correlated to the instance size as it is the case for the random instances. Since some of the instances here exhibit several primary funnels, note that f go represents the relative size of largest of those funnels. Other instance features seem to influence the landscape global structure. The number of funnels on the structured instances is somewhere in between that of the uniform and the Table 3: Random instances global structure metrics. Relative sizes of the primary funnel f go , the largest funnel f lg , proportion of solutions in more than one funnel f ov , number of attractor local optima a, number of funnel sinks s, and relative in-strength of global optimum sinks d go . Chained-LK success rates are reproduced in the last column for comparison purposes.  clustered random instances (see Table 3). The uniform random instances show the largest number of funnels in the studied set. For the instances with high levels of neutrality, i.e. nrw1379 and the drilling problems except u574 and rat783, it was not possible to assess whether the attractor local optima at a given objective function level fully connected into a single plateau or groups of plateaus sharing the same objective function value. The plateaus are potentially very large, with the sampling process not guaranteeing their full exploration. Alternatively there may be multiple plateaus for a single objective function level. We therefore used two different methods to estimate bounds for the different metrics. The first method approximates the plateaus by assuming that they are indeed con-nected. With this approach, instances such as d1291 and fl1577, reveal a relatively small number of funnels despite their large number of different local optima. The second method does not assume that there is a single plateau for each objective function level, instead it considers all connected nodes sharing the same evaluation as different plateaus. This leads to 1000 such plateaus, the same as the number of runs, for instances u1060 and fl1577. Further discussion of this issue, through the analysis of distances between solutions, is presented in Section 5.4. Note that for these two instances, f go and f lg have very small ranges and are thus summarised by a single value.
A thorough study of these highly neutral instances may require additional sampling efforts and additional analysis to characterise the extent of the plateaus. We leave that for future work. Nevertheless, the global structure in terms of the number of funnels does not seem to differ significantly from that of the city instances or clustered random instances, if we assume that plateaus are indeed connected. We argue that search difficulty on very neutral landscapes relates not only to the multifunnel structure, but also to the size and location of the plateaus. A large plateau at the global optimum level may reflect an easy to search instance, while a large suboptimal plateau may act as a trap that is difficult to escape from.

Network visualisation
One of the advantages of modelling a system as a network, is the possibility of visualising it. This is one of the strengths of the proposed approach, as it allows a more accessible way of grasping the complexities of landscapes global structure.
Software for analysing and visualising networks is currently available in various languages and environments. Here we use the R statistical language together with the igraph package (Csardi and Nepusz 2006). Layout algorithms are at the core of network visualisation, they assign vertices to positions in a metric space. Force-directed methods model the pairwise attraction and repulsion of vertices, and are known to reflect the community structure or modularity of a network (Noack 2009). We, therefore, use them in order to visually characterise the landscapes' multi-funnel structure. Funnels can be visually identified as modules in the network. As our model indicates, nodes are LK-search local optima and edges represent escape transitions according to double-bridge moves. We decorated them to reflect features relevant to search dynamic. The size of nodes is proportional to their incoming strength (weighted incoming degree), therefore, it reflects the extent to which nodes attract the search dynamics. The colour of nodes reflects their funnel membership. We used the heat colours palette, a sequential colour scheme skewed to the reds and yellows. Red identifies the global optima, and the yellow colour gradient reflects an increase in cost. The edges' widths are proportional to their weight, which indicates the frequency of transitions. That is, the most frequently visited edges are thicker. We present both 2D and 3D images. In Ochoa and Veerapen (2016a), we proposed a 3D visualisation where the x and y coordinates are, as usual, determined by a graph layout algorithm; the innovation is to use the objective function as the z coordinate. This provides a clearer representation of the funnel sink and basin concepts, bringing an almost tangible as-pect to the landscape and funnel metaphors. The global optimum can be identified in the 3D plots as the node with the lowest z coordinate.
In order to have manageable images, we plotted the networks corresponding to the subset of local optima within 0.1% or 0.05% in evaluation from the global optimum. We also removed self-loops for improved visibility. Figures 4 and 5 illustrate local optima networks for selected random and structured instances, respectively.
The plots in Figures 4a and 4b illustrate how the number of funnels rapidly increases with the number of cities for the uniform random instances. Instance E1243.85 features 76 and 19 funnels for solutions within 0.1 and 0.05% in evaluation from the global optimum, while instance E755.73, have 4 and 2, respectively. The size of the funnel containing the global optimum is much smaller for E1243.85 as compared with the smaller uniform instance E755.73 and the clustered instance of the same size C1243.85. Indeed, instances E755.73 (Fig. 4a) and C1243.85 (Fig. 4b) have a similar Chained-LK success rate despite their difference in size, which indicates the importance of global structure in search difficulty. Both instances feature a large secondary funnel coloured in orange, which acts as a strong attractor of the search process as indicated by the size of its funnel sink. Figure 5 illustrates the local optima networks for solutions within 0.1% in objective function value from the optimal solution for 2 city instances att532 and pr1002, and one drilling problem u574. Instances att532 (Fig. 5a) and u574 (Fig. 5c) have a similar CLK success rate and also a similar global structure for solutions within 0.1% in evaluation from the global optimum cost. They both feature two funnels visualised in red (the primary funnel) and yellow (the secondary funnel). In att532 the two funnels are overlapping, whereas in u574 they are separated (for solutions within 0.1 % in evaluation, but they may overlap at higher solution costs). An important difference between the random (Fig. 4) and the structured instances (Fig. 5) is that the latter show neutrality. This is reflected by the number of global optima (a single one for all the studied random instances), and 2 and 4 for att532 and u574, respectively. The neutrality can also be appreciated at higher solution costs on the 3D plots for att534 and u574, where several solutions are located at the same objective function level.
The global structure of city instance pr1002 (Fig. 5b) is strikingly different. It clearly has a large primary funnel sink visualised in red, which reflects the high CLK success rate (0.67) of this instance. There are however 3 secondary funnels whose sink solutions have an evaluation 80 units higher than the optimum (visualised in orange and yellow). It is important to note that these 3 nodes do not form a plateau despite having the same objective function values, they are not connected with doublebridge moves according to our sample. These 3 funnels have several connections to other solutions in the primary funnel, although they do not connect directly to the primary funnel sink. Therefore, they are separate funnels according to our empirical definition. The visual depiction however, seems to reflect that they belong to the primary funnel, which suggests that there might be hierarchies of funnel structures.
For some of the studied instances, each sampling run visits a different set of solutions. This leads to as many connected components as there are runs, as is the case for instances u1060 and fl1577. Furthermore, each component consists of a long chain of nodes with some small loops when plateaus are explored. The previously used network visualisation provides relatively little information for such cases.  presents an alternative visualisation for instances where different runs do not share any common solutions. Blue dots represent single solutions. Red lines represent consecutive solutions in a plateau, i.e. the solutions have the same cost, and the length of the line is proportional to the number of solutions in the plateau. Each run is displayed as one column of dots and lines since there is no overlap between the solutions found in different runs. Figure 6 specifically shows the top 25 largest components, i.e. runs, containing the 25 largest plateaus. A marked difference can be observed in the distribution of solutions and plateaus between the two instances. For u1060, there is usually a single plateau at the end of the run that is very close to the global optimum. The latter is not actually reached for the top 25 runs with largest plateaus. For fl1577, runs usually encounter a number of different plateaus, manage to escape from them, and finally get stuck in a non-optimal plateau. The globally optimal value is reached in only 3 out of 25 runs. The fact that u1060 has fewer plateaus to escape than fl1577 provides an explanation for its higher success rate in finding a global optimum out of 1000 runs. The visualisation provides a straight-forward view of the distribution of plateaus.
A question that remains is whether all the plateaus with equal evaluation found across the different runs are actually a single very large plateau. We attempt to answer this by analysing the distance between solutions in Section 5.4.

Distance analysis
In addition to visualising and analysing the networks, it is also useful to examine the landscapes using the pairwise distances between nodes. We consider the bond distance in particular, which is defined as the difference in the number of common edges, or bonds, between two solutions. It corresponds to the number of edges within a solution minus the number of edges that are in common between the two solutions considered.
Bond distance heatmaps can be useful to quickly assess how close pairs of solutions are to each other. Figure 7 displays the pairwise bond distance between sink nodes or nodes within sink-plateaus for the six instances visualised in Figures 4 and 5. The objective function of each node is displayed on the vertical axis. The plot is mirrored along the diagonal. For E1243.85 and pr1002, only the ten best local optima are considered and, in both cases, these nodes are sink nodes and are not part of a plateau. For the other four instances, there are ten or fewer nodes to consider. The colour gradient allows us to clearly distinguish the two plateaus within instances att532 and u574. The bond distance within these plateaus is lower than 15 units and generally only in single digits, showing that the solutions are very similar, as expected.
When considering highly neutral instances, like u1060 and fl1577, where each run produced a new connected component, an important question is whether the different plateaus found at the same objective function value are actually part of some very large plateau or if they are an artefact of the stochasticity of the sampling method. Figure 8 attempts to provide some insight by visualising the distribution of pairwise bond distances between sink-plateau nodes sharing the same evaluation. In practice, because each sink-plateau contains several hundred nodes, we choose the node   Fig. 7: Heatmaps for pairwise bond distance between sink nodes or nodes within sink-plateaus. Labels on the y-axis indicate objective function but omitted on the x-axis. Only the 10 best nodes are considered for E1243.85 and pr1002. For both, the nodes considered are sink nodes and are not part of a plateau. The two sinkplateaus in att532 and u574 are clearly identified by the nearly-white shading between neighbouring nodes (low bond distance).
with the highest density as the representative of a given sink-plateau. Ties are broken randomly. Each representative node is then compared to the other representative nodes that share the same evaluation. The plot has bond distance as x-axis, frequency of occurrence as y-axis and the objective value of each node pair is colour-coded as shown in the legend. For the two instances, the distributions appear multimodal with one major peak. This was confirmed using Hartigan's Dip test (Hartigan and Hartigan 1985;Maechler 2016) which rejected the null hypothesis of unimodality (p-value < 2.2 × 10 −16 ). The smallest bond distance found for u1060 is 17, while it is 42 for fl1577. The largest bond distances are 193 and 199 for u1060 and fl1577, respectively. The distances are fairly large, meaning that it is not easy to move between plateaus sharing the same evaluation. Still, for u1060, the smaller distances may indicate that these is a roundabout way to bridge the disconnected plateaus. For fl1577, the larger distances indicate it is unlikely that the presence of multiple sink-plateaus is purely the result of the sampling process. However, we cannot conclusively confirm that the sinkplateaus found are actually separate either. There might be some pathways between sink-plateaus with the same evaluation if plateaus were explored for more than 10 000 iterations.

Effects of sampling methodology
In this subsection we investigate the influence of the stopping criterion and of the number of runs on the networks and their metrics. In a sense, these parameters affect the "depth" and "width" of the sampled networks.
We first consider the impact of using fewer than 10 000 consecutive iterations without improvement as the stopping criterion, and whether this threshold is sufficient for the different metrics to converge. For 18 of the 20 instances, we generate networks for 1000 to 10 000 iterations without improvement, with steps of 1000. The instances left out are those with high neutrality, u1060 and fl1577, for which 1000 runs are not enough to find common solutions between runs. For these instances, a new sampling methodology that considers the extensive neutrality needs to be devised, but this is left as future work.
For each (metric, instance) couple, values outside the [0, 1] range are normalised and, for each stalling iteration increment, the average is computed across the instances. Results are presented in Figure 9 where each line joins the average values across instances for each metric. The success rate is highlighted by the red line. Some of the curves are markedly non-monotonic. This is an artefact of using independent subsets of runs for each stalling threshold. In addition, the normalisation exacerbates the differences between those subsets.
As can be observed, on average the different metrics have approximately converged when 10 000 non-improving iterations are used as stopping criterion. For many metrics the majority of the change in values occurs before the 5000 iterations mark. For a few (a, a p and s p ), this happens later, around the 9000 iterations mark. The two weighted degree metrics (d and d max ) naturally do not converge since more connections between already discovered nodes appear, especially in sink plateaus. Overall, these findings suggest that the stopping criterion is appropriate and the sampled landscapes would not differ greatly with a higher threshold. We repeated a similar experiment for the number of runs with values from 100 to 1000 with steps of 100. The stopping criterion is fixed to 1000. Normalisation was performed where necessary as described earlier. We discarded a few data points for which no global optimum was found since several of our metrics take the global optima into account. Results are presented in Figure 10.
In this scenario, the metrics fall within two broad categories: ones that are relatively constant and those whose value increases. The latter are metrics that depend on the size of the sample: for instance, the number of local optima and unique evaluations, the in-strength, and the maximum length of a path to a global optimum. The number of clusters and of attractor and sink plateaus also increases on average across the instances. However, when looking at instances individually, they have converged by 1000 runs for smaller instances such as att532. Let us note that the number of global optima remains more or less constant, except for those instances that have a lot of them.
Metrics that maintain relatively constant values are those describing properties of the sampled network with respect to the whole: the proportion of nodes for which there is a path to a global optimum, the average length of such a path, the relative sizes of funnels and the relative in-strength of global optimum sinks. On the whole, these metrics might seem preferable since they do not depend on the number of runs. Nevertheless, metrics where convergence with respect to the number of runs is unlikely may still be useful and some arbitrary number of runs should be chosen, probably related to the available computing budget.

Correlation study
In order to obtain more general insights into how different landscape features affect search difficulty, we conduct a correlation study of the different metrics studied in  this paper with respect to the success rate across 1000 runs. To assess the relative quality of local optima network and global structure metrics as estimators of success, we also include 64 features based on TSP instance characteristics. These include features that describe the edge cost distribution, cluster characteristics and minimum spanning trees (MSTs). They are computed using the tspmeta R package (Mersmann et al 2013). Pearson correlations are calculated across two separate sets of instances: 100 uniform and 100 clustered instances. For the sake of clarity, only the top 5 positively and negatively correlated features are presented in Table 5.
We found both instance and network structure features that strongly correlate with Chained-LK success rate. The features showing the strongest correlations differ between both sets of instances. On clustered instances, a higher number of structural features show a strong correlation while the opposite is observed on uniform instances. However, on the two instance sets, the average path length to a global optimum,l go , is strongly correlated to success. Also, d go , the normalised incoming strength of global-optimum sinks and its complement, d ngo , the incoming strength of non-globally optimum sinks show strong correlations. This is because sub-optimal sinks act as traps to the search process, from where Chained-LK search cannot escape with its perturbation operator. For clustered instances, several of the top correlations are related to neutrality (n, evals, n/evals).
When considering all network structure metrics discussed previously in this paper, almost all of them produced correlation values below −0.4 or above 0.4. The features with low correlation values for clustered instances are the number of connected components (c), the proportion of solutions in more than one funnel ( f ov ) and the mean and maximum in-strength (d and d max ). In contrast, for uniform instances, the d max shows a strong correlation (0.74) and it is only f ov that exhibits a low correlation value with success. Overall, this suggests that the metrics studied reflect the search dynamics of Chained-LK.

Conclusions
Revealing what makes a combinatorial problem hard remains an open challenge. In this quest, understanding the global structure of the underlying landscapes is essential. There are few attempts in the literature to analyse, let alone visualise, the global structure of combinatorial landscapes. This is due in part to the lack of adequate tools to study their complexity. Local optima networks help to fill this gap. In our sample of TSP problem instances, we found evidence of multiple funnels, instead of a single big-valley as previously believed, in TSP landscapes of moderate size (500 to over 1500 cities). Good local optima decompose into multiple valleys of different depths, each channelling the search process to a separate low cost solution or group of solutions. We also found evidence of high amounts of neutrality or extensive plateaus at the local optima network level in several of the structured instances (i.e. city and drilling problems).
Our data-driven approach models realistic landscapes as networks, empirically characterises the notion of funnels according to notions from theoretical chemistry, and proposes novel 2D and 3D network visualisations. This brings new quantitative and visual insights into landscapes' global structure. The 3D plots provide a concrete and intuitive illustration of the fitness landscape and funnel metaphors. The depiction of edges incoming strength as node sizes, visually reveals the strong attractors of the search process. We also propose alternative visualisation tools for analysing very neutral, and conduct a distance analysis. Finally, we conduct a correlation study between Chained-LK success rate and both landscape and TSP instance features.
We found significant differences among the studied instance classes. Randomly generated instances have a single global optimum and reveal null or very small neutrality. On the other hand, structured instances from the TSPLIB generally have more than one different global optima, with some instances featuring a large global optimum plateau. This is probably because there are many pairwise city distances with the same value. Within the same instance class, the number of funnels generally increases, while the size of the funnel containing the global optima generally de-creases with the size of the instance. However, the instance class strongly influences the global structure. The random uniform instances revealed a much larger number of funnels as compared with the other instances studied, which explains why these instances are generally harder to solve with heuristic search methods.
Our correlation study reveals that Chained-LK success is strongly correlated to several landscape structural features. On clustered instances a higher number of structural features show high correlations, as compared to the uniform instances where a large number of instance features show strong correlations. For both instance types, however, the features showing the strongest correlation are the incoming weighed degree of the global optimal sink and its complement, the incoming weighed degree of the sub-optimal sinks. This confirms the importance of landscapes global structure in explaining search difficulty: sub-optimal sinks act as traps to the search process, from where Chained-LK cannot escape with its perturbation operator.
The impact of neutrality on search difficulty is harder to assess, and this will motivate future work with additional sampling mechanisms to explore the extent of local optima plateaus. Preliminary experiments adding a stronger perturbation to Chained-LK proved to help in smoothing the funnel structure, that is, reducing the number of funnels and making the global optima more reachable (Ochoa and Veerapen 2016a). Future work will expand this study, explore the role of crossover operators in landscapes with multiple funnels, and apply the methodology to other combinatorial optimisation problems.