{'fraction_non_alphanumeric': 0.03774353191755936, 'fraction_numerical': 0.00967863183612103, 'mean_word_length': 4.20662100456621, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fast Algorithms for Indices of Nested Split Graphs Approximating Real Complex Networks

We present a method based on simulated annealing to obtain a nested split graph that approximates a real complex graph. This is used to compute a number of graph indices using very efficient algorithms that we develop, leveraging the geometrical properties of nested split graphs. Practical results are given for six graphs from such diverse areas as social networks, communication networks, word associations, and molecular chemistry. We present a critical analysis of the appropriate perturbation schemes that search the whole space of nested split graphs and the distance functions that gauge the dissimilarity between two graphs.


Introduction
The architecture of real-life networks may not be predictable and involves data which is sometimes difficult to manage. Algorithms known to date for the computation of many parameters associated with complex networks are often exponential-time, blocking further development. To overcome this problem, the strategy implemented here is to use the class of nested split graphs (NSG) since a graph in this class can provide a model network which is sufficiently close to a given real-world network and which enables easy computation. The class of NSGs provides an ideal candidate as these graphs are nice to work with, are easily stored, and their structured topology lends itself to polynomial time algorithms for the computation of certain of its invariants. Where this is not possible, their rich mathematical versatility may enable other methods to circumvent the problem.
Nested split graphs, also known as threshold graphs, form a subclass of split graphs in which the vertex set is partitioned into a clique (maximal complete subgraph) and a co-clique, that is an independent subset of vertices (with no edges between any pair of the independent vertices). Furthermore, NSGs display a structure in which the clique and the co-clique are partitioned into cells which form an equitable vertex partition, that is each vertex of a cell has the same number of neighbours in each of the cells. Moreover, the adjacencies among the cells impose a structure on the NSG that induces a partition of the vertices referred to as the NSG partition.
The questions we ask are: (a) given a network G, is there a NSGĜ that is sufficiently close to G that can be used instead of G for the purposes of computing network parameters and spectral properties? and (b) how can computation time be reduced by drawing on the geometric properties of NSGs?
To search for the optimumĜ, we devise two distance functions, the scaled walks and the spectral functions, f S W and f λ , respectively, with the objective to measure the rate of convergence of the original given network to an interim NSG G picked by the algorithm considered. The techniques used are Markov chains and simulated annealing that enable the space of coded NSGs to be searched. Three perturbation schemes are used for each of the distance functions. These are the 'hamming', 'edge', and 'move' schemes, all of which will be explained in detail in Section 3.2. The process outputs the best network among all NSGs that can replace the given network, up to a prescribed tolerance, for the purpose of computing selected invariants associated with the original network.
The parameters to be computed using the different techniques are derived mainly from sociology, information theory and physical chemistry. They are Entropy, the Randić, Wiener, Szeged, Co-PI and Estrada indices, and also Gutman's Graph Energy and Resolvent Energy. The networks tested come from communication networks, word associations, and molecular chemistry.
When using Markov chains and simulated annealing to determine a NSG close enough to a given G, the problem centres on the best distance function and Markov chain perturbation to use for reasonable computation times and reliable values of graph invariants of the NSG to which the process converges. The different techniques are analysed to see which distance functions and perturbations correspond to optimized output and well behaved algorithms.
The paper will be organised as follows. In Section 2 we start by presenting various characterization of NSGs. This is followed in Section 3 by descriptions of the graph parameters to be estimated, and the distance functions to be used as a measure of closeness between the interim NSG G in the simulated annealing process and the original graph G. The simulated annealing process is presented in Section 4, together with the perturbation schemes used to traverse the whole space of NSGs. In this section we also develop efficient algorithms to determine the graph parameters of Section 3 by taking advantage of the geometric properties of NSGs. Estimates of selected parameters of G are obtained by calculating them forĜ, the limiting NSG. Practical results are given in Section 5 for a selection of six graphs, where we consider the effect of the simulated annealing parameters on the NSGs obtained, and also compare the indices for the limiting NSGs with those of the original graphs. In Section 6 pointers for further work from these seminal ideas are suggested.

Characterizations of a Nested Split Graph
A complete graph on n vertices is K n and a graph on n vertices with no edges is K n . The path and the cycle on n vertices are denoted by P n and C n , respectively. A disconnected graph consisting of r copies of a graph G is denoted by rG. The number of edges incident to a vertex v of a graph G on n vertices is the degree ρ v . A graph has an isolated vertex if it has a vertex with no edges incident to it, while a dominating vertex corresponds to a vertex with edges to all the other vertices of the graph. A vertex v is a duplicate of a vertex u in a graph G if u and v are not adjacent and they have the same neighbouring vertices in G. Two vertices u and v are co-duplicates if they are adjacent and have the same neighbours.
Different characterizations of NSGs which we shall use in computing the different NSG parameters are now presented.

Theorem 1
The following statements are equivalent for a graph G [1,2]: (ii) G is P 4 , C 4 and 2K 2 free; (iii) G is constructed from K 1 by successive additions of an isolated or dominating vertex.
From Theorem 1(iii), the following result gives the construction of a connected NSG from a binary code.
Proposition 1 A connected NSG on n vertices is uniquely coded as a binary string of (n − 2) bits. The rth bit is 0 if the (r + 1)th vertex added is an isolated vertex and it is 1 if it is a dominating vertex. This string, which is referred to as the minimum representation of the NSG, assumes that the first vertex in the construction is an isolated vertex and that the last vertex is a dominating vertex.
Observe that (n − 2) bits are sufficient because in the canonical construction the first vertex is always isolated and the last vertex has to be dominating. It also follows from Theorem 1(iii) that each intermediate subgraph in the construction of the NSG G is itself a NSG. If the first bit of the encoding string (reading from the left) is 0, then the intermediate NSG G 2 on two vertices is 2K 1 , whereas if it is 1, then G 2 is K 2 .  It is convenient to make use of a more concise form of the creation sequence by encoding repetitions of successive bits.

Definition 2
The compact creation sequence for a connected NSG on n vertices is the r-tuple a = (a 1 , a 2 , . . . , a r ), where a i = n, r is even, and for all i, a i ≥ 1 is the length of the ith successive sequence of entries in the creation sequence of the same value (i.e. with each entry equal to 0 or each entry equal to 1).
Recall that a clique or a co-clique in the NSG vertex partition is referred to as a cell. Note that the compact creation sequence is a direct representation of the NSG partition with the smallest even number r of cells. Therefore, the coefficients of a with even subscripts give the size of the cliques, while those with odd subscripts give the size of the co-cliques in the NSG. By considering the construction of the n-vertex NSG from the compact creation sequence a = (a 1 , a 2 , . . . , a r ), the next result follows.
Proposition 2 If to a NSG G on n vertices, a duplicate to a vertex in the co-clique or a co-duplicate to a vertex of the clique is added, then a NSG G on n + 1 vertices is obtained with the same number of cells.
The adjacency matrix A = (A i, j ) of a NSG can be obtained directly from the creation sequence c as Observe that there is no computation involved in this process.

Measures
Given a graph G = (V, E), we would like to find a nested split graph (NSG)Ĝ that is sufficiently close to G. Closeness in this sense implies thatĜ has similar properties to G where it matters within a given context. Such properties are usually defined in terms of one or more indices of the graph. We review the indices on which the algorithms are based.

Indices of a graph
We consider first indices based on graph topology. The simplest of these is a function of the number of edges, |E|. Specifically, for an undirected graph G = (V, E) where all edges have unit weight, the entropy [3][4][5][6] is given by Increasing in complexity, we consider indices that are a function of the vertex degrees. The Randić index [7] of a graph G = (V, E) is given by where ρ u and ρ v are the degrees of vertices u and v respectively. A further increase in complexity is necessary for indices that are a function of the shortest distance between vertex pairs. The simplest of these is the Wiener index [8-10] of a graph G = (V, E) given by where d(u, v) is the shortest distance between vertices u and v. A related metric is the Szeged index [11,12] of a graph }| is the number of vertices that are closer to u than to v. A variation on the Szeged index, the Co-PI index [13] of a graph G = (V, E) is given by Next we consider indices based on the graph's spectral properties. The Estrada index [14,15] of a graph G = (V, E) is defined by where n = |V| and λ j , j = 1, . . . , n are the eigenvalues of the graph's adjacency matrix A. A related metric is Gutman's graph energy [16][17][18][19], obtained as which is the sum of absolute values of the eigenvalues of the graph's adjacency matrix. Finally, the resolvent energy [20,21] is also obtained from the eigenvalues of the graph's adjacency matrix, and is defined as

Measures of distance between graphs
The problem of finding a NSGĜ that is sufficiently close to a given graph G can be seen as a combinatorial optimization problem, where the objective function to be minimized is a measure of the distance between the two graphs. This raises the inevitable question of what objective function should be used. One approach is to choose a graph parameter and use either the absolute difference or the Euclidean distance between the parameters of G and G. An alternative is to choose a vector function that represents the properties of interest for a graph, and use say the Euclidean distance between the vectors representing G andĜ as a measure of dissimilarity. We propose two vector functions, representing the graph's walk and spectral properties.
The first vector function W is based on the number of walks from each vertex, and represents the reachability of vertices from each other. Specifically, for a graph G = (V, E) with adjacency matrix A we define the scaled walk matrix by which is the concatenation of scaled walk vectors W i for walks of length i, given by where ∆ = max v∈V ρ v is the maximum vertex degree of G and j is the all-one vector of length n = |V|. The second vector function Λ represents the graph's spectral properties. Specifically, for a graph G = (V, E) with adjacency matrix A we use the vector of eigenvalues [22][23][24] defined by where n = |V| and λ j , j = 1, . . . , n are the eigenvalues of A indexed in non-increasing order.

Simulated annealing
Simulated annealing is a probabilistic global optimization algorithm that mimics the physical process of thermal annealing, by which crystals are slowly cooled to reach a state with the lowest free energy. In analogy, the algorithm seeks the lowest energy state from a large but finite state space, where the energy is a deterministic function of the state.
The algorithm traverses a Markov chain, where the probability of accepting a state change depends on the difference in energy ∆ E between the current and candidate states, and on a time-varying parameter T called the temperature. Specifically, when the candidate state has a lower energy than that of the current state, then ∆ E < 0 and the acceptance probability is 1, while when the candidate state has an energy equal to or higher than that of the current state, this is accepted with probability e − ∆ E T . The initial temperature needs to be large enough so that a transition to a higher-energy state is allowed with high probability. As the temperature is reduced, such transitions are accepted with lower probability. The final temperature is chosen so that only transitions that reduce the energy are allowed with any significant probability. Typically the annealing schedule follows a geometric decrease from the initial to final temperature.

Perturbation schemes
Recall that for a NSG, the creation sequence c = (c 1 , c 2 , . . . , c n ) must satisfy the constraints c 1 = 0 and c n = 1. Therefore, we define the state variable that uniquely corresponds to this NSG as the tuple σ = (c 2 , . . . , c n−1 ).
The transition probabilities of the state space define a neighbourhood, where ideally the neighbours of a state have an energy similar to that of the initial state. The neighbourhood must be defined in such a way that, while traversing the Markov chain, all states are reachable, following a sufficient number of transitions. This ensures that the whole state space can be explored. The neighbourhood can be defined as a perturbation function on the state variable. We propose three different perturbation schemes: 'hamming', 'edge', and 'move', that define neighbourhoods within the state space.
With the 'hamming' scheme, the neighbourhood is defined as the set of those nested split graphs with a valid creation sequence at a Hamming distance of one from the current state's creation sequence. The neighbouring states are those corresponding to creation sequences c = (c 1 , c 2 , . . . , c n ), where and j ∈ {2, . . . , n − 1} identifies the specific neighbour. Observe that this definition creates a path between all possible states, with a minimum distance of at most n − 2 steps. In our implementation, the neighbourhood is uniformly sampled.
Observe also that with the 'hamming' scheme not all neighbouring graphs are equally similar to the current state. For example, consider the current state corresponding to c = (0, 1, 1, 1, 1), equivalent to K 5 . One possible neighbour corresponds to c = (0, 0, 1, 1, 1), which has one less edge than c, while another corresponds to c = (0, 1, 0, 1, 1), which has two fewer edges than c. In order to minimise the differences between neighbours, we define the 'edge' scheme as the perturbation whose neighbourhoods form a subset of the 'hamming' neighbourhood where the changed coefficient is the first or last variable coefficient or it is different from either of its adjacent coefficients. Thus The latter condition allows the movement of a vertex from a cell to another cell with adjacent subscript, that is from a clique 2k to a co-clique 2k − 1 or 2k + 1 and from a co-clique 2k − 1 to a clique 2k or 2k − 2. This includes cases where, as a result of the change, the number of vertices in the initial cell is reduced to zero, and forces the consequent merging of both cells that are adjacent in index to the vanishing one. This condition provides paths to all states with equal or smaller number of cells than the current state. The former condition allows the inversion of coefficients c 2 and c n−1 , and provides paths to states with an increased number of cells.
Finally, we define the 'move' scheme, where neighbouring states are obtained from the current one by moving one of the vertices in a source cell j to a different target cell k where k is at most two different from j. Formally, for a current state defined by the compact creation sequence a = (a 1 , a 2 , . . . , a r ), neighbouring states are those corresponding to compact creation sequences a = (a 1 , a 2 , . . . , a r ), where and 1 ≤ j ≤ r, j ∈ N is the cell a vertex is moved from while 1 ≤ k ≤ r, k ∈ { j − 2, j − 1, j + 1, j + 2} is the target cell. This condition provides paths to all states with equal or smaller number of cells than the current state. The neighbourhood also includes states with r + 2 cells, obtained by moving two vertices from cell j to form two new cells of size 1. These states correspond to compact creation sequences a = (1, 1, a 1 , a 2 , . . . , a r ) and a = (a 1 , a 2 , . . . , a r−1 , a r , 1, 1), where and a j ≥ 2, j ∈ {1, 2} in the case of a , while j ∈ {r − 1, r} in the case of a . Note that the above definition includes degenerate compact creation sequences, where the size of one cell becomes zero. In these cases, both cells that are adjacent in index to the vanishing one will merge.

Fast computation of indices on NSGs
The structure of NSGs lends itself to less complex algorithms for computing a number of graph indices. We consider here the indices introduced in Section 3.1, describing algorithms to compute these for NSGs, and comparing the complexity of these algorithms with those for general graphs.
We start with the computation of entropy. This requires the computation of the number of edges in a graph, which for a general simple undirected graph can be computed as: where A i, j is the element at row i, column j of A. This is a summation over n 2 −n 2 elements, for a complexity O(n 2 ). For a NSG with compact creation sequence a = (a 1 , a 2 , . . . , a r ), the number of edges can be computed as a summation of two terms using the elements of a: where κ k = r/2 j=k a 2 j . The first term gives the number of edges in the subgraph containing all cliques, while the second term considers the edges between each co-clique and all the cliques it is connected to. This requires the computation of the cumulative sum of even-indexed elements followed by a summation over the odd-indexed elements of a, for a complexity O(r). The computation of entropy follows directly, with constant complexity.
Consider next the Randić index, which is a function of the degrees of the vertices at either end of an edge. Observe that all vertices within a cell have the same degree. For co-clique a 2k−1 on level k ∈ {1, . . . , r/2} the degree is equal to the sum of vertices in cliques at the same level or higher: Similarly, for clique a 2k , the degree is equal to the sum of vertices in co-cliques at the same level or lower plus the number of other vertices in cliques: Each vertex degree requires computation of complexity O(r), for a total complexity O(r 2 ) to compute all the vertex degrees. The Randić index can then be obtained by summing over the edges between cells and the edges within cliques, taking into account the multiplicity in each case: Respectively, the terms count over the edges between co-cliques and cliques, edges between distinct cliques, and edges within cliques. Computationally, this can also be improved to where α k = a k √ ρ k and β k = r/2 j=k α 2 j . This involves four computations of complexity O(r), reducing the complexity from O(r 2 ) when (21) is used.
The Wiener index introduces a further increase in complexity, requiring the shortest distance between every vertex pair. For a general undirected graph, this can be obtained with a breadth-first search starting from every vertex, for an overall complexity O(n 3 ). For a NSG we can observe that: (a) the distance between two vertices in cliques is always equal to one, (b) the distance between two vertices in co-cliques is always equal to two (they are always connected to at least one common clique vertex), (c) the distance between a vertex in a co-clique and a vertex in a clique at the same level or higher is always equal to one, and (d) the distance between a vertex in a co-clique and a vertex in a clique at a lower level is always equal to two. Now the Wiener index can be computed as the sum of the shortest distance between each distinct pair of vertices, which leads to As before, this can be simplified to where γ k = r/2 j=k a 2 j−1 and δ k = r/2 j=k a 2 j . Again this reduces the computation to complexity O(r). The Szeged index first requires the computation of ν(u, v), the number of vertices that are closer to u than v, for every directly connected pair of vertices u, v in the graph. Now for a NSG there are three types of directly connected vertices: (a) when both in the same clique, (b) when the vertices are in different cliques, and (c) when one vertex is in a co-clique and the other is in a clique at the same level or higher. Consider these in turn, observing that ν(u, v) ≥ 1 as there will always be at least one vertex (u itself) that is closer to u than v. If u, v are in the same clique, they are co-duplicates; so it follows that ν(u, v) = ν(v, u) = 1. Without loss of generality, if u is in clique 2k and v is in clique 2 j where k < j and k, j ∈ {1, . . . , r/2}, then u has the same connections as v except for co-cliques a 2i−1 , i ∈ {k + 1, . . . , j}, where v has a direct connection while u does not. It follows, in this case, that ν(u, v) = 1 and ν(v, u) = 1 + j i=k+1 a 2i−1 . Finally, if u is in co-clique 2k − 1 and v is in clique 2 j where k ≤ j and k, j ∈ {1, . . . , r/2}, then we need to consider the distances from u, v to other vertices w in the graph: (i) if w is in a co-clique at the level of v or lower (including the same co-clique as u), it is adjacent to v but not to u, so that d(v, w) < d(u, w), except when w = u, in which case d(u, w) < d(v, w) (ii) if w is in a co-clique at a level above v, it is adjacent to neither u nor v, and d(v, w) = d(u, w) = 2, (iii) if w is in a clique at the level of u or higher (including the same clique as v), it is adjacent to both u and v, so that d(v, w) = d(u, w) = 1, except when w = v, in which case d(v, w) < d(u, w), and (iv) if w is in a clique at a level below u, it is adjacent to v but not to u, so that d(v, w) < d(u, w). Combining these for case (c), it follows that ν(u, v) = 1 and Finally, to obtain the Szeged index we need to sum over all edges of the graph: For computational purposes, this can be simplified to where k = k j=1 a 2 j−1 , ζ k = k j=1 a 2 j , η k = r/2 j=k a 2 j j , and δ k is as defined earlier. This involves a sequence of computations of complexity O(r).
The Co-PI index is also based on ν(u, v); following earlier observations for the Szeged index, we can similarly obtain the Co-PI index as a sum over all edges in the graph: As for the Szeged index, this can be simplified to where k , ζ k , η k , and δ k are as defined earlier. Again, this involves a sequence of computations of complexity O(r).
The remaining indices are functions of the eigenvalues of the graph's adjacency matrix A. These eigenvalues are referred to as main if they have some eigenvector not orthogonal to the all-one vector j, and non-main otherwise [1]. Now for a NSG it can be shown that the main eigenvalues of A can be obtained by computing the eigenvalues of Q where Q is the r × r matrix whose rows are the distinct rows of AX, for X being an n × r indicator matrix satisfying AX = XQ [1]. The matrix Q can be obtained directly from the compact creation sequence a as Unlike A, Q is not a symmetric matrix, so the computation of its eigenvalues is more complex by some constant factor. However, this is more than made up for by the significant reduction in size from n × n to r × r. This reduces the complexity of computing the eigenvalues from O(n 3 ) to O(r 3 ). In addition to the eigenvalues of Q, A has two distinct non-main eigenvalues: (i) the eigenvalue −1 for co-duplicates, with multiplicity r/2 k=1 (a 2k − 1) = δ 1 − r 2 , and (ii) the eigenvalue zero for duplicates, with multiplicity r/2 k=1 (a 2k−1 − 1) = γ 1 − r 2 . We can therefore obtain the Estrada index as a summation whereλ j , j = 1, . . . , r, are the eigenvalues of Q. Similarly, Gutman's graph energy is obtained as and the resolvent energy is obtained as

Graphs used
Results in this section are for a number of simple connected graphs. We start with two small graphs taken from social communities, the results for which can be more easily verified. These are named G12 and G14, respectively with 12 and 14 vertices, 13 and 18 edges, and topologies as shown in Figures 2 and 3 respectively. Next we consider the graph 'geant', representing the GÉANT network topology as of January 2017 [25]. This has 41 vertices and 64 edges, with topology as shown in Figure 4. We also consider the graphs 'abstract' and 'honey', representing word associations in text paragraphs. These have 59 and 77 vertices, 605 and 687 edges, respectively. Finally we consider the graph 'benzenoid', representing the carbon bonds of the large benzenoid known as the supernaphthalene made by Klaus Müllen [26]. This has 96 vertices and 129 edges, with topology as shown in Figure 5.

Simulated annealing
For each of the graphs described in Section 5.1, we use simulated annealing to find a similar nested split graph, by minimising one of the distance functions. Two graph distance functions are used: 1. the walk matrix distance, which is the Euclidean distance between the two graphs' scaled walk matrices, and 2. the spectral distance, which is the Euclidean distance between the two graphs' eigenvalues. For each case, we also use three different perturbation schemes as defined in Section 4.2. Throughout we use an annealing schedule where the temperature is varied logarithmically from 10 2 to 10 −7 over one million steps.    The final output of the simulated annealing algorithm is tabulated in Table 1 for graphs G12, G14, 'geant', 'abstract', 'honey', and 'benzenoid'. Comparing the behaviour of the graph distance functions as optimization targets, we note that the spectral distance is more stable, often converging to the same graph for different perturbations. The most significant exception is for the graph 'honey'; perhaps this is because 'honey' is not connected while all the other graphs are connected. On the other hand, the walk matrix distance converges to different graphs, except for the smallest cases (G12 and G14).
When we consider the effect of the perturbation scheme on the simulated annealing algorithm, we can draw the following conclusions. The 'hamming' scheme consistently performs best, resulting in the lowest energy case every time there is a difference between methods, with the exception of 'benzenoid' where the convergent distance for the three methods is very similar. This comes at a cost: since the 'hamming' neighbourhood contains a greater number of NSGs that are significantly different from the current state, the acceptance rate during annealing falls more rapidly. In turn, this requires a slower reduction in temperature, as shown later in the annealing timeline. As expected, 'move' is consistently better than 'edge', as it has a substantially wider neighbourhood, without including substantially different NSGs. We can also see that 'edge' favours more compact creation sequences, while 'hamming' results in longer ones. This can be explained by observing that neighbours in 'edge' generally have the same number of cliques and co-cliques. An increase only happens with an insertion at either end of the creation sequence, while a decrease only happens with a move from a clique or co-clique with one vertex. As expected, the performance of 'move' is somewhere in between.
It is also helpful to compare the annealing timeline under different conditions. These are shown in Figure 6 for the graph 'geant'. For the range of temperatures considered, annealing starts with a 100% acceptance rate and 50% improvement rate, indicating that the algorithm is correctly exploring neighbourhoods indiscriminately. The acceptance rate starts to drop as the temperature is decreased, with a corresponding drop in improvement rate, until both reach zero. At lower temperatures, the algorithm explores less and favours only changes towards the local minimum. This behaviour is reflected in the energy values sampled at different temperatures, which exhibit large variation at higher temperatures and converge towards the final low energy state at lower temperatures. Timelines for the other graphs exhibit a similar pattern.

Comparison of indices for similar graphs
Finally we compare the values of a number of indices on graphs for the best output using each graph distance function with those of the original graph. These are tabulated in Table 2 for graphs G12, G14, 'geant', 'abstract', and 'benzenoid'. Note that we do not include 'honey' because the original graph is not connected, so that many of these indices cannot be computed. From these results we can observe the following general patterns. Limiting NSGs optimized for spectral distance tend to have a number of edges that is closer to that of the original graph than those optimized for walk distance. Consequently, indices that directly depend on edge count (such as entropy) are better approximated by these graphs. These graphs are also better approximations for indices that depend directly on the eigenvalues, as is to be expected. On the other hand, graphs optimized for walk distance are better approximations when computing indices that depend on the distribution of vertex degrees, such as the Randić index. In fact, the approximation here is often particularly close. Finally, for indices that depend on the distribution of distances between vertex pairs, the better approximator depends on the specific index. The Wiener index and Co-PI index are better approximated by graphs optimized for spectral distance, while those optimised for walk distance are better for the Szeged index.

Conclusions
In this paper we have presented a method based on simulated annealing to obtain a NSG that approximates a given real complex graph. Three different perturbations have been considered for traversing the space of NSGs on a given number of vertices, and two distance measures to determine how well a specific NSG approximates the original graph. Practical results for six graphs from social networks, communication networks, word associations, and molecular chemistry, show that the simulated annealing process is stable and converges to suitable NSGs.
We have also developed very efficient algorithms to compute a number of useful graph indices, taking advantage of the geometrical properties of NSGs. For the same six graphs, we have compared the indices of the original graphs with those of the corresponding approximate NSGs, using each of the two distance functions. We see that the ideal distance function to use depends on the indices to be computed.
The results shown in this paper are for relatively small graphs, with at most 100 vertices and 1000 edges. For practical purposes it would be interesting to see how well the methods presented here scale to much larger graphs, such as those from the Stanford Large Network Dataset, with 10-10 000× more vertices and edges. The efficient algorithms presented here to compute the graph indices scale well to graphs of that size, because their complexity is polynomial in the number of cells of the NSG, rather than polynomial in the number of vertices. This is a consequence of the structural properties of the NSG. However, the simulated annealing method as presented here does not scale so well. Optimizing the simulated annealing method is possible, and we leave this for future work.