On Near Optimal Spectral Expander Graphs of Fixed Size

09/29/2021 ∙ by Clark Alexander, et al. ∙ 0

We present a pair of heuristic algorithms. The first is to generate a random regular graph of fixed size. The second is the introduction of the Metropolis Coupled Simulated Annealer (MCSA) for optimizing spectral gaps in fixed size regular graphs.



There are no comments yet.


page 6

page 14

page 15

page 18

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

From the 1980s and 1990s one has seen several methods for generating random regular graphs [MW,W]. While generating a random graph is extremely easy, the constraint of requiring regularity is what produces the difficulty. There is, however, some renewed interest in random regular graphs as we see that many random regular graphs become near spectral expanders [BB]. In addition, regular graphs allow for consistent edge labeling in a rotation map. [A2] which allows one to efficiently compute a discrete time quantum walk [A1].

This work is presented in two major parts. The first is presenting an extremely simple method for producing a random regular graph with a fixed number of vertices and degree of regularity. This method involves a constructive method for building a regular graph of fixed size and a randomization algorithm by which we can move around a fixed number of edges.

The second part which entails the lion’s share of the present work shows how to produce a Metropolis Coupled Simulated Annealing algorithm for near global optimization of spectral gaps in random regular graphs of fixed size. In the initial stages of this work, we had tried to use a simple simulated annealer, which in large part was able to find Ramanujan graphs quickly and easily. However, during our computational explorations, we found that the method we used for picking a “neighboring” solution greatly affected the overall quality of the solution. Furthermore, switching the method of choosing a neighbor during the annealing process again affected the quality of solution. Thus the next level of complexity called for coupling the annealing solutions. This draws us back to a Metropolis Coupled Markov Chain Monte Carlo (MCMCMC) [MNL] which entered into the computational space cerca 1997 as a way of having a more effective Monte Carlo method for multi-modal distributions. We have modified the MCMCMC to be a Metropolis coupled simulated annealer (MCSA) to look for near optimal spectral expander graphs. We give some numerical results and some empirical evidence in support of an open question as to whether there are infinitely many 7-regular Ramanujan graphs. In short, we have tested every graph with on even number of vertices from 8 to 1000 vertices, and have always been able to produce a Ramanujan graph within a matter of seconds. We have also tested graphs of size 1000 to 2000 adding 20 vertices at a time and have again been able to find Ramanujan graphs, but the computational effort grows substantially after approximately 1200 vertices. Nonetheless, we have found several 7-regular 2000 vertex Ramanujan graphs and the adjacency matrix for one such is available upon request from the author.

2 Part 1: The Algorithms for Randomization

2.1 Constructing a Regular Graph of Fixed Size

In order to draw a regular graph we only need to check one criterion, that is, the product of number of vertices with degree of regularity is even. In our case we will only consider simple undirected graphs with no loops. From our first week of graph theory class we remember for a graph with degree of regularity we have.


We cannot satisfy this if and

are both odd.

Thus we can split our regular graph generation into two cases for the parity of . In the first, and easier case, we let be even. This allows us to pick any number of vertices (greater than ). The first step is to label our vertices to ( or 0 to mod ). And create the cycle graph . This gives us a 2 regular graph. Now we iterate for to connecting to . Below is pseudo which has been implemented in both Python(3) and Julia (1.4). This is a pseudo-code for generating the adjacency matrix of a regular graph. We assume that the product of vertices and degree of regularity is even.

1:procedure GenerateRegularGraph(vertices, degree))
2:      = vertices
3:      = degree
4:     if  then
5:          This corresponds to connecting diametrically opposed vertices.
6:         for  in  do
9:         end for
10:          Symmetrize the matrix
11:     else We don’t have diametrically connected vertices
13:         for  in  do
16:         end for
18:     end if
19:     return
20:end procedure
Algorithm 1 Rotation map from the adjacency matrix; Pseudo-code

Let’s take a look at three potential examples. These are all produced by the networkx package in Python. Figure 1 gives 4 regular graphs with no randomness the number of vertices and degrees of regularity are respectively.

(a) 9 vertex, 4-regular
(b) 12 vertex 4-regular
(c) 12 vertex 6-regular
(d) 14 vertex 5-regular
Figure 1: Regular graphs with no randomness.

2.2 The Algorithm for Randomization

The general algorithm is quite straight forward once we have a regular graph in place. In plain English we algorithm looks something like:

  1. Pick two edges at random, swap one end point of each.

  2. If the graph is still regular with the same degree and there are no loops accept the switch.

  3. If not, pick again.

  4. Repeat until “sufficiently” randomized.

The supporting procedures of checking for loops (i.e. the diagonal of the adjacency matrixhas a nonzero element) and that the graph is still regular (all rows sum to the same degree) will be assumed. Now we come to the pseudo-code implementation of the graph randomization algorithm. First, we need the ability to switch edges. In the adjacency matrix representation the edges are represented by locations in the matrix which are 1. So the procedure getEdges is simply find the nonzero locations in the adjacency matrix

1:procedure SwitchEdges(vertices, degree)
2:     (vertices, degree)
3:      In python we need .copy() to get a different memory pointer
4:     edges getEdges()
5:     Randomly choose two edges from edges.
16:     if is_regular() and np_loops(then
17:         return
18:     else
19:         return
20:     end if
21:end procedure
Algorithm 2 Switch Edges; Pseudo-code
1:procedure GenerateRandomGraph(vertices, degree, switches)
2:     (vertices, degree)
3:      SwitchEdges
4:     for  in 1:switches do
5:          SwitchEdges
6:     end for
7:     return
8:end procedure
Algorithm 3 Generate Random Regular graph; Pseudo-code

2.3 Examples of Randomized Regular Graphs

Now we give three pairs of examples. In order to show that we can use this to generate random graphs of nearly arbitrary degree. We’ll give examples of size (20,8), (30,13), and the much larger (500,25).

(a) 20 vertex 8-regular nonrandom
(b) 20 vertex 8-regular random
Figure 2: 20 vertex 8-regular graphs
(a) 30 vertex 13-regular nonrandom
(b) 30 vertex 13-regular random
Figure 3: 30 vertex 13-regular graphs

For the larger example, the plots are not so enligthening, but we shall show them here and give timing so that we can demonstrate the overall efficacy of our ability to randomize graphs of nearly arbitrary degree.

(a) 500 vertex 25-regular nonrandom
(b) 500 vertex 25-regular random
Figure 4: 500 vertex 25-regular graphs

The second graph in 4 is the results of 1500 random edge switches and on the author’s laptop took 53.6s to complete.

And for posterity let’s examine the spectral propertires of these two graphs.

(a) 500 vertex 25-regular nonrandom
(b) 500 vertex 25-regular random
Figure 5:

500 vertex 25-regular graphs eigenvalue histograms

It’s clear from the spectrum that both graphs are 25 regular and that the randomized graph is closer to the semi-circular distribution than the nonrandom graph.

3 Part 2: The Search for Ramanujan Graphs which do not have Explicit Constructions

3.1 How to Anneal a Regular Graph

By the construction in §2.1 we achieve a regular graph of fixed size. In particular we build the adjacency matrix of a regular graph of fixed size. This construction, however, has a small spectral gap and a large graph diameter. Using our randomization scheme we wish to construct an adjacency matrix of a regular graph with a larger spectral gap. We have seen before [BB] that a random regular graph comes close to a Ramanujan bound for spectral expansion. The goal at present is to use our randomization scheme to search for ever larger spectral gaps. In reality, we will use a simulated annealer to minimize the normalized . Looking a little further down 8 shows how eigenvalues perform with single random switches.

We see that that randomization scheme first brings to within the Ramanujan threshold quickly, but eventually moves above it. This does not, however, get to the weak bound of Alon-Bopanna, nor the strong bound. It is our goal to see how close we can get to these bounds. The striking feature or many graphs is how quickly we can cross the simple Ramanujan threshold.

The general structure of a simulated annealing algorithm is straight-forward:

  1. Guess a solution at random

  2. Score that solution

  3. Pick a “neighbor” solution

  4. Score the neighbor.

  5. If the neighbor is better, make it the new solution

  6. if neighbor is worse, still accept neighbor if not too far off

  7. Lower the algorithm’s “temperature”

  8. Repeat until satisfied

This is a gross oversimplification, but there are two salient points we need to consider. The first being: How do we pick an effective “neighbor” solution. The second being: How do we deal with the temperature parameter? Changes in either of this points can drastically affect the quality of solutions in a simulated annealer. In particular, lowering the temperature too quickly (i.e. quenching) tends to land a solution in the first deep local minimum the algorithm finds. On the other hand, lowering the temperature too slowly, takes too long and begins to look like a brute-force solution, which we don’t have time for.

In the case of the random regular graph, we will set our annealing algorithm to look for the lowest possible normalized rather than the largest spectral gap. As it is set up, however, we will never achieve a bipartite graph as the second normalized eigenvalue is 1. In particular a bipartite graph of degree has eigenvalues

, which we normalize and consider absolute values. Thus this algorithm ignores bipartite graphs as good expanders since we wish to achieve the steady state vector quickly by starting a random walk at a single vertex rather than in a superposition of vertices in different partition sets. We see, from the figure

8 that switching a pair of edges does well, however, for larger and sparser graphs, this takes too long. Consider, for example, a large cubic or 4-regular graph. That is a graph with several thousand vertices. A single edge switch does very little to affect the eigenvalue gap. Thus, one may consider using multiple edge swaps as the neighboring solution. This means that a “neighbor” graph should swap a dozen or so pairs of edges before measuring .

As with a standard annealing algorithm we want to move from essentially a random walk when the temperature is hot to a stochastic gradient descent type algorithm. We, however, wish to slow down not only the random selection of temporary solutions, but also slow down or speed up the edge switching as necessary. Thus we have decided to move from a simple single simulated annealer to a coupled annealer.

3.2 Metropolis Coupled Simulated Annealing

A single simulated annealer is an example of a Markov chain Monte Carlo method tuned for optimization as opposed to finding a posterior probability distribution. As discussed in the previous section we wish to couple multiple simulated annealing algorithms to find a solution which is closerto a global maximum. The idea of coupling was first (published) in the late 1990s [MNL] with an eye toward bimodal and multimodal distributions where MCMC algorithms are slow to converge. From the perspective of this work, we shall leverage multiple optimizers to look for a (closer to) global optimum. In order to tune our coupling we implement the following procedure:

  1. Initialize simulated annealers with initial temperatures with different cooling rates These are called chains.

  2. Each chain gets a different method of selecting a neighboring solution. For this work the chain has a neighbor selection of random edge switches.

  3. For each chain run one step of the simulated annealing algorithm and cool each chain

  4. Rank order the chains by best solution. The best current solution gets the lowest temperature. This is the “cold chain.”

  5. Randomly select two other chains, score them and compute an acceptance probability for switching (or more simply one can simply exchange them without the extra computation).

  6. Select a number from a uniform variable . If the acceptance probabilty is higher than the randomly drawn number, switch the chains. That is, switch their neighbor solutions, temperatures and cooling rates.

  7. Repeat steps 3-6 until the hottest chain has reached the temperature threshold, or we have reached an acceptable bound (e.g. We have found a Ramanujan graph, or reached a weak lower bound for ).

3.3 Some Bounds and Asymptotics

For this section we will present the results as they were originally stated and then show the results for normalized eigenvalues. We will also let a graph have degree of regularity throughout.

3.3.1 Eigenvalues

The original bound for an expander graph to be a Ramanujan graph is

In our case, of course, this means

Additionally from [N, W1] we have several other bounds. There is a weak bound in which we have


where . This particular is chosen so that where

The basic construction comes from finding an eigenvector corresponding to a particular vertex within

. With a small inspection one can see that as this bound becomes the original Ramanujan bound. In our case, we will look for eigenvalues which cross the (slightly) higher bound

We will demonstrate that for smaller graphs(up to roughly 2000 vertices) we can often get below this bound.

There is one additional much stronger bound


where is the diameter of the graph. In [N] this is a strict lower bound. Again, one can see as we again approach the Ramanujan bound. This bound as we will show empirically can not be crossed, and even approaching this bound is difficult except for complete graphs or particularly tiny graphs. For example, we have enough computing fire power that we can check every regular graph of size .


We can precompute the value of the diameter () of the graph by a simple counting argument. This is done by assuming we have a -ary tree. The distance from the root to any leaf is the radius which is half the diameter.


3.3.2 The Number of Regular Graphs

For smaller graphs we have exact calculations of the number of regular graphs [OEIS1, OEIS2]. We also have in [MW] (p58 theorem 1) the number of -regular graphs on vertices is


Rearranging this a bit we get the logarithmic count of -regular graphs on vertices is


3.3.3 Cycles in Random Regular Graphs

We see in [BB, W2] that the number of cycles of length in a

-regular random graph is asymptotically a Poisson random variable

Pois with mean


3.4 Known Constructions

There are several well-known constructions for families of Ramanujan graphs, in particular, LPS, Morgenstern, Paley, MSS, Zig-Zag, and Super-Singular Isogeny graphs. Additionally [F,BB] points out that with high probability a -regular random graph is Ramanujan or near Ramanujan (ie )

The known constructions cover the following cases

Construction Degree of Regularity Conditions
Margulis-Gabber-Galil [Mar, GG] 4 vertices or infinite
Morgenstern [Mor] for prime
Zig-Zag [RVW] are degrees of regularity of the two graphs to be multiplied
Paley [Pa, ER] or prime, prime
SuperSingular Isogeny Graphs [Pi, EHLMP] prime.
Bipartite Graphs [MPS] any Graph must be Bipartite

We see now that the smallest positive integer for which a family of constructions is not known is for graphs with chromatic number 3 or more. We also have a few other small unknown cases Additionally, we see that there are some graphs with specific numbers of vertices which are not known in constructions. For example, the Zig-Zag product requires the number of vertices to be composite and Paley graphs have vertices which are elements of a finite field so we have powers of primes. In the case where the smallest degree of regularity for a Paley graph is which quickly becomes quite large. We will explore some relatively smaller graphs for visualization purposes and some moderate sized graphs for the purposes of numerical demonstration.


The methods of [Mar,GG] cover 4-regular graphs, but with a perfect square number of vertices. The MCSA finds a 17 vertex Ramanujan graph which is 4-regular in using a flag to stop at finding such a graph. Additionally, the MCSA finds a near optimal graph of 17 vertices () in . The Ramanujan threshold for a 4-regular graph is . Among 17 vertex 4-regular graphs there are according to [MW]. This is a small demonstration of the efficacy of MCSA. In the next section we shall explore some numerical results more deeply.

3.5 Some Numerical Results

To demonstrate the efficacy of the heuristic search involved, let’s begin with some well-known graphs which are near optimal, the fullerenes [AKS]. In particular we’ll look explicitly at the dodecahedron and the “bucky-ball.” These are well known 3-regular (cubic) graphs with 20 and 60 vertices respectively. To get our eyes adjusted to the way we will draw graphs, here are depictions of the respective fullerene graphs.

Figure 6: The Dodecahedron and the Bucky Ball

In our scenario we initialize two graphs, a 20 vertex 3-regular and a 60 vertex 3-regular which look like this.

Figure 7: 3-regular graphs from explicit construction

These graphs are nowhere near expanders, much less optimal expanders. In particular the normalized threshold for a cubic graph to Ramanujan is . In our case the initialized graphs have for the 20 vertex graph and for the 60 vertex graph.

Again, from [F,BB] we know with high probability that a random regular graph will be Ramanujan. So let’s follow two progressions randomly switching edges for each of these graphs and see how often they fall below the required threshold.

Figure 8: Eigenvalues of 20 and 60 vertex graphs via random switching

A couple of immediate observations. The smaller the graph, the easier it is to be a spectral expander. This agrees with our intuition about edge expanders as we require very few moves with a small number of vertices to reach any vertex via random walk when the graph is “small.” Larger graphs, however, require more moves and we see this in the fact that a single random switch affects far less in the 60 vertex case than in the 20 vertex case. We also see that a vast majority of random regular graphs are below the Ramanujan threshold. In particular in two random trials of 499 graphs, we have 345 such graphs below the threshold in the 20 vertex case and 432 in the 60 vertex case. Additionally we have 68 and 196 graphs below the bucky ball line respectively. However, nothing gets below the dodecahedron line.

Consider now two graphs that we have generated with our Metropolis coupled simulated annealer. These are noth first run algorithms, and thus if we had tuned our parameters very specifically, we could have achieved lower eigenvalues.

Figure 9: Expanders of 20 and 60 vertex graphs via MCSA

Again we notice that neither of these graphs can quite hit the dodecahedon threshold, but they have both surpassed the bucky ball threshold easily as well as exceeding any graph achieved by single random switches.

Now for just a moment let’s turn our attention to 7-regular graphs. This is one of the cases which is not explicitly covered in known constructions. Drawing a graph with more than 60 or so vertices becomes rather unattractive on a page this small, so we will stick to 50 vertices for the sake of demonstration. In a timed run, we flagged the algorithm to stop once it had a achieved a Ramanujan graph. This particular 50 vertex, 7-regular graph cost the author’s laptop 1.1373 seconds of cpu time. The normalized threshold for a 7-regular graph is

Figure 10: A quick Ramanujan graph with

As we can see, this is significantly lower than the required Ramanujan threshold. This is even below a weak optimal expander as defined in Equation 2.

Now consider a graph that we did not flag.

Figure 11: A near optimal graph with

In the case of 7-regular graphs the thresholds we’re dealing with are

weak optimal
strict lower bound

In nearly every case attempted by the author the heuristic optimizer was able to find a graph close to (and in many cases below) the weak optimal expansion constant. In every case, we were able to find a Ramanujan graph, but again, this is not surprising based on the work of [F] However, we were never able to find a graph which met the strict lower bound as given by [N]. In simple numerical trials we were able to find Ramanujan graphs with degree of regularity 7 for every (even) number of vertices between 8 and 2000. This work, however, does not answer the question of whether there is an explicit construction for such graphs.

While the drawing of a graph with 2000 vertices is not very attractive, we can still show the histogram of normalized eigenvalues to show the spectral expansion. Below is a histogram of a 2000 vertex 7-regular graph discovered by Metropolis coupled simulated annealing. The author has also added a spreadsheet of its adjacency matrix into github (email author for address) so that researchers can explore its random properties further at their desire.

Figure 12: Histogram of 2000 vertex 7-regular Ramanujan graph normalized eigenvalues;

Appendix: Additional Algorithms

1:procedure PartialAnneal(matrix, trials, temperature, cooling_rate, switches)
2:      = matrix
3:      = temperature
4:      = computeSecondEigenvalue
5:     solution0 = A
6:     best_ =
7:     best_solution = solution0
8:     for k in 1:trials do
9:         solution1 = getNeighbor(solution0, switches)
10:          = computeSecondEigenvalue(solution1)
11:         if  best_ then
12:              best_ =
13:              best_solution = solution1 Keep the global best in storage
14:         end if
15:         if  then
16:              solution0 = solution1
17:               =
18:         else
19:              ap = acceptance_probability(, , )
20:         end if
21:         if ap rand() then
22:              solution0 = solution1
24:         end if
25:     end for
26:     cooling_rate
27:     return best_solution, best_,
28:end procedure
Algorithm 4 Partial Annealing
1:procedure DefineCoupling((vertices, degree, chains, min_cooling, max_cooling))
2:     graphs = [generateRegularGraph(vertices, degree) for k in range(chains)]
3:     initial_temperatures = ones(chains)
4:     cooling_rates = [start = min_cooling, end = max_cooling, spacing = (max_cooling - min_cooling) / chains]
5:     return graphs, initial_temperatures, cooling_rates
6:end procedure
Algorithm 5 Define Coupling
1:procedure PerformOneStep(graphs, temperatures, cooling_rates)
2:     _list =
3:     for k in 1: length(graphs) do
4:         graphs
5:          = temperatures
6:         _decay = cooling_rates
7:         , = PartialAnneal(, , cooling_rate = _decay, switches )
8:         _list _list
9:         temperatures =
10:         graphs
11:     end for
12:     ReverseRankOrder = argsort_list
13:     graphs = graphs[ReverseRankOrder]
14:     temperatures = temperatures[ReverseRankOrder] Now randomly swap a pair of graphs
15:     if length(graphs)  then
16:         randomLocations = sampleFrom([1:length(graphs)], number = 2, without replacement)
17:         swap(graphs[randomLocations[1]],graphs[randomLocations[2]])
18:     end if
19:     return graphs, temperatures, minimum(_list)
20:end procedure
Algorithm 6 Perform One Step
1:procedure CoupledAnnealing(vertices, degree, chains, min_cooling, max_cooling, _min)
2:     graphs, s, cooling_rates = DefineCoupling(vertices, degree, chains, min_cooling, max_cooling)
3:     best_ = computeSecondEigenvalue(graphs[1])
4:     best_ = graphs[1]
5:     while maximum(s) _min do
6:         graphs, s, current_ = performOneStep(graphs, s, cooling_rates)
7:         if current_ best_ then
8:              best_ current_
9:              best_ = graphs[1]
10:         end if
11:     end while
12:     return best_, best_
13:end procedure
Algorithm 7 Coupled Annealing