Would a grid-like graph have lower average shortest distances than a random graph?

Question

Imagine two different types of street networks with the same number of nodes and where edges are weighted according to their length. Both network types cover the same geographic area (say 1 km²) and have the same node density. The first type is the grid, where every street goes in one of two possible, perpendicular directions and where every street is the same length (a). The other type corresponds to more random networks, where streets don't follow a particular order, their length and their directions varying randomly (b and c, for instance). In each network, you have eight points of departure/arrival, labelled A to H. You measure the average shortest path distance from every point to every other point. Would the grid network have shorter shortest paths on average?

Intuitively, I would assume not. I would say that the grid network would have less variable shortest paths, since you don't risk a big detour, but you are also unlikely to find a shortcut. However, I was looking to see if there are any graph theoretical/topological reasons why a regular network might be more efficient than a more random one.

Answer 1

Summary

I wrote a Python script to numerically search for the optimal solutions of your problem (see below). While this may not find the optimum, it does find solutions that are better than the grid. On the other hand, an arbitrary random network within any reasonable confines is very unlikely to beat the grid.

How the script works and why it is not perfect

The script considers the average shortest path as a function of the position of the freely movable nodes. It then numerically tries to find a minimum of that function using SciPy’s local and global optimisation routines.

We here want to do a global optimisation in a high-dimensional parameter space ($2n$ dimensions, where $n$ is the number of free nodes). The problem here is that there are lots of local minima and we do not know where to search. Also, to evaluate a given parameter set, we have to compute the corresponding shortest paths, which takes time. We therefore need to keep the number of these evaluations small.

You can roughly compare this to trying to find the deepest point of the oceans with lead and line. It is easy to find the local minimum, i.e., the lowest point of the valley you are in: You just measure depth again a few hundred metres on either side and go for the steepest descent (that is roughly what SciPy’s minimize does). However, to find the global minimum, you first need to find the trenches to begin with, and that requires you mapping the entire ocean first. And if you can do so only coarsely, and do not have any knowledge of tectonic plates, you may easily miss the Mariana trench and just find some other rather deep point.

Your example (2×2)

The optimal solution for 2×2 free nodes is almost certainly this one (or a rotation thereof):

with an average shortest path length of 3.076, while it is 3.142 for the grid. (No, that’s not π, but $\frac{22}{7}$.)

Bigger networks

For bigger networks, it’s even more difficult to be sure about the absolute optimum, but it’s never the grid. Below are the best solutions I could find for the respective sizes. Interestingly, they become less and less symmetric.

Script

import numpy as np
from matplotlib.pyplot import subplots
from scipy.optimize import minimize, basinhopping, dual_annealing
from scipy.sparse.csgraph import csgraph_from_dense, floyd_warshall

class Network(object):
    def __init__(self,shape):
        self.shape = np.asarray(shape)
        self.build_topology()
        self.prepare_graph_and_parsing()
        # initialise with a grid:
        self.apply_parameters(np.ones((*self.shape,2)))
        # array containing the relevant entries of the path matrix for quicker averaging:
        self.relevant_paths = np.tril_indices(len(self.fixed_nodes),k=-1)

    def build_topology(self):
        # nodes are identified with their positions in a grid
        edges = set()
        self.free_nodes = list()
        self.fixed_nodes = list()

        for i in (0,self.shape[0]+1):
            for j in range(1,self.shape[1]+1):
                self.fixed_nodes.append((i,j))
        for i in range(1,self.shape[0]+1):
            for j in (0,self.shape[1]+1):
                self.fixed_nodes.append((i,j))

        for i in range(1,self.shape[0]+1):
            for j in range(1,self.shape[1]+1):
                self.free_nodes.append((i,j))
                for neighbour in [ (i+1,j), (i-1,j), (i,j+1), (i,j-1) ]:
                    edges.add( ( (i,j), neighbour ) )

        self.nodes = self.fixed_nodes+self.free_nodes
        self.edges = list(edges)

    def prepare_graph_and_parsing(self):
        # initialise locations as perfect grid, so the fixed ones stay like this
        self.locations = np.indices(self.shape+2,dtype=float).transpose(1,2,0)

        n = len(self.nodes)
        node_to_idx = { node:i for i,node in enumerate(self.nodes) }
        weight_matrix = np.full((n,n),np.nan)
        for k,(A,B) in enumerate(self.edges):
            # fill weights with edge number so we can sort by it
            weight_matrix[node_to_idx[A],node_to_idx[B]] = k
        self.graph = csgraph_from_dense(weight_matrix,null_value=np.nan)

        # Sort edges by appearance in graph.data to make things easier
        self.edges = [ self.edges[int(k)] for k in self.graph.data ]

        # Create arrays indicing the components of the edges, so they can be easily extracted from the location array:
        edge_pos_gen = lambda: np.empty( len(self.edges), dtype=int )
        self.edge_1_x = edge_pos_gen()
        self.edge_1_y = edge_pos_gen()
        self.edge_2_x = edge_pos_gen()
        self.edge_2_y = edge_pos_gen()
        for k,(A,B) in enumerate(self.edges):
            self.edge_1_x[k],self.edge_1_y[k] = A
            self.edge_2_x[k],self.edge_2_y[k] = B

    def apply_parameters(self,parameters):
        # `parameters` are the advances from the previous point in respective row/column
        # This avoids intersections.
        self.parameters = parameters.reshape((*self.shape,2))
        for i in range(2):
            self.locations[1:-1,1:-1,i] = np.cumsum(self.parameters[:,:,i],axis=i)

        # Calculate the distances for all edges and directly feed the graph.
        self.graph.data = np.linalg.norm(
                  self.locations[self.edge_1_x,self.edge_1_y]
                - self.locations[self.edge_2_x,self.edge_2_y]
            , axis=1)

    def avg_path(self):
        paths = floyd_warshall(self.graph,directed=False)
        return np.average(paths[self.relevant_paths])

    def badness(self,parameters):
        self.apply_parameters(parameters)
        return self.avg_path()

    def optimize(self,**kwargs):
        return minimize( self.badness, x0=self.parameters, **kwargs )

    def optimize_bh(self,**kwargs):
        self.apply_parameters( basinhopping(
                self.badness,
                x0 = np.ones(self.parameters.size),
                **kwargs,
            ).x )

    def optimize_da(self,**kwargs):
        n = self.parameters.size
        self.apply_parameters( dual_annealing(
                self.badness,
                bounds = np.vstack(( np.full(n,0), np.full(n,3) )).T,
                x0 = np.ones(n),
                **kwargs
            ).x )

    def visualise(self,filename):
        fig,axes = subplots()
        axes.set_xlim(0,self.shape[0]+1)
        axes.set_ylim(0,self.shape[1]+1)
        axes.set_xticks(range(self.shape[0]+2))
        axes.set_yticks(range(self.shape[1]+2))
        axes.set_aspect("equal")
        axes.set_title(f"average shortest path: {self.avg_path():.5f}")

        for X,Y in self.edges:
            axes.plot(*zip(self.locations[X],self.locations[Y]),"blue")

        for X in self.free_nodes:
            pos = self.locations[X]
            axes.annotate( f"$\\binom{{{pos[0]:.3f}}}{{{pos[1]:.3f}}}$", pos, color="red" )

        fig.savefig(filename,bbox_inches="tight")

minimiser_kwargs = {"method":"Powell", "options":{"ftol":1e-5}}

for shape,iterations in [ ((2,2),1000), ((3,2),1000), ((3,3),1000), ((4,2),1000), ((4,3),300), ((4,4),100), ((5,5),10) ]:
    print(shape)
    my_network = Network(shape)
    my_network.visualise(f"initial_{shape[0]}_{shape[1]}.png")

    my_network.optimize( **minimiser_kwargs )
    my_network.visualise(f"solution_{shape[0]}_{shape[1]}.png")

    my_network.optimize_bh( T=1, minimizer_kwargs=minimiser_kwargs, niter=iterations, stepsize=1 )
    my_network.visualise(f"solution_{shape[0]}_{shape[1]}_bh.png")

    my_network.optimize_da( local_search_options=minimiser_kwargs, maxiter=iterations, initial_temp=5e4, visit=2.99 )
    my_network.visualise(f"solution_{shape[0]}_{shape[1]}_da.png")