Strongly regular graphs (SRG), Paley, Brouwer-Haemers graphs.

[KeesP]

## brouwer-haemers-graph.r
## Unlicense (Ø) 2023 CLP, MR-Amsterdam.
## -------------------------------------
## Strongly regular graphs SRG(v, k, λ, μ).
## Wikipedia (en), Strongly regular graph.
## Paley graphs,          SRG(v = 4t + 1, k = 2t, λ = t − 1, μ = t). 
## Brouwer Haemers graph. SRG(v = 81, k = 22, λ = 1, μ = 6).  
##
## References.
## Andries E. Brouwer & Hendrik Van Maldeghem, Strongly regular graphs, Amsterdam and Ghent, March 2021.
## Ch. 7.4.4 Paley graphs, Ch.10.28 The Brouwer-Haemers graph.
## [🔗](https://homepages.cwi.nl/~aeb/math/srg/rk3/srgw.pdf)
## 
## Brouwer, A.E.; Haemers, W.H, Structure and uniqueness of the (81, 20, 1, 6) strongly regular graph, Tilburg University.
## [🔗](https://pure.uvt.nl/ws/portalfiles/portal/996688/structur.pdf)
##
##  Brouwer, A.E., Paley graphs.
## [🔗](https://www.win.tue.nl/~aeb/graphs/Paley.html)
##
## Brouwer, A.E., Brouwer-Haemers graph
## [🔗](https://www.win.tue.nl/~aeb/graphs/Brouwer-Haemers.html)
##
##  Code layout.
## (i)   Theory, modern algebra
## (ii)  Finite field operations and structures
## (iii) Creating finite fields
## (iv)  Calculating residuals
## (v)   Creating and plotting the graph
## (vi)  Statistics
## (vii) Examples   
## 
## ----------------------------------------------------------------- (i)
## Some modern algebra (Artin, Galois Theory, p28).
## q is prime. Zq, or Z/qZ, integers modulo q, %% q in R, Zq is a field (K) iff q is prime.
## https://math.stackexchange.com/questions/252046/the-paley-graph-with-9-vertices
##
## q is a prime power p^d. GF(m)[X]/P is a field iff P is a prime ideal.
## Let K[x] be the polynomials with coefficients in K.
## Let p(x) be polynomial anx^n + an-1x^n-1 ... a0.
## Let p(x)K[x] be ideal P, consisting of all polynomials with factor p(x).
## Let F = K[x]/P be the ring of cosets. A field iff p is irreducible in K.
## Let q be the order of F. Note that P is the zero-element in F.
## F is homeomorphic to K(ξ) the field of formal expressions an-1ξ^n-1 ... a0.ξ^0.
## Powers of ξ less then n are linearly independent and
## ξ^n ≡ -(an-1ξ^n-1 ... a0.ξ^0).
##  
## An element is represented in R as vector c(an-1, ... a0).
## ξ^k corresponds with position q-k in vector c(...).
## Multiplication is implemented as shift register (to the left).
## Let p(x) = x^2 + x + 2, coefficients in K = Z3, -1 ≡ 2.
## ξ^2 ≡ -1ξ - 2, or  ξ^2 ≡ 2ξ + 1 this generates the cyclic field F9:
## ξ^0 = 1, ξ1 = ξ, ξ2 = 2ξ + 1, ξ3 = 2ξ + 2, ξ4 = 2, ξ5 = 2 ξ6 = ξ, ξ7 = ξ + 2, ξ8 = ξ + 1.
##
## Irreducible, primitive polynomials over GF(3)
## https://baylor-ir.tdl.org/bitstream/handle/2104/8793/GF3%20Polynomials.pdf?sequence=1
##
## ------------------------------------------------------------------- (ii)
## Primitive polynomials over finite fields GF(q).
## p(x) = m               GF(m), m is prime, d=1.
## px  <- c()
##
## Rook's graph, generelized quadrangle, GQ(2,1).
## p(x) = x^2 + x + 2     GF(3^2). m=3, d=2.
## px  <- -c(1, 2)        # primitive polynomial in K.
##
## p(x) = x^2 + x + 2     GF(5^2). m=5.
## px  <- -c(1, 2)        # primitive polynomial in K.
##
## p(x) = x^2 + x + 2     GF(7^2). m=7.
## px  <- -c(1, 3)        # primitive polynomial in K.
##
## p(x) = x^2 + x + 7     GF(11^2). m=11.
## px  <- -c(1,7)         # primitive polynomial in K.
##
## p(x) = x^2 + x + 2     GF(13^2). m=13.
## px  <- -c(1,2)         # primitive polynomial in K.
##
## p(x) = x^3 + x + 2     GF(5^3). m=5, d = 3.
## px  <- -c(1,1,3)       # primitive polynomial in K.
##
## BrHa                   Brouwers Haemers graph when TRUE
## p(x) = x^4 + x + 2     GF(3^4). m=3, d=4.
## px  <- -c(0, 0, 1, 2)  # primitive polynomial in K.
## -----------------------------------------------------------------------
BrHa  <- TRUE             # Adjacent vertices differ nonzero fourth power.
m     <- 3                # K = GF(m), E = GF(m^d).
px    <- -c(0,0,1,2)      # Primitive polynomial in K. 
d     <- length(px)       # d is degree of extension.

## calculate powers of ξ^n, ξ^n+1
pwx <- list()
nxt <- rep(0, d); nxt[1] <- 1     
for (i in seq_len(d-1) ) {
   nxt      <- (c(tail(nxt, d-1), 0) + head(nxt, 1) * px) %% m
   pwx[[i]] <- nxt
}

## multiply by ξ and reduce
gf_shift <- function(x, px, q=m) { n <- length(px); (c(tail(x, n-1), 0) + head(x, 1) * px) %% q }

## multiply in quotient ring K[x]/P)
## polynom package
## https://stackoverflow.com/questions/39979884/calculate-the-product-of-two-polynomial-in-r
gf_mult <- function(a, b, px=px, q=m) {
  n  <- length(a)
  ab <- outer(a, b)
  p  <- as.vector(tapply(ab, row(ab) + col(ab), sum))

  ## reduce powers x^q, x^q+1
  for (i in seq_len(n-1)) {
    p[- (1:n-1)] <-  tail(p, n) + pwx[[i]] * p[n - i]
  }
  return(tail(p, n) %% q)
}

## ---------------------------------------- (iii)
## create Galois field of order q (GFq).
## points in GF(q) == graph vertices.
if (d > 1) {
  start <- rep(0, d); start[d] <- 1
  alpha <- c(rep(0, d-2), c(1,0))
  nxt   <- start

  GFq <- list(rep(0, d))                # Galois field including zero.
  for (i in seq_len(m^d)+1) {
    GFq[i] <- list(nxt)
    nxt <- gf_mult(nxt, alpha)          # Shift left and reduce, multiply by ξ.
    if (all(nxt == start)) break        # Full cycle iff polynomial is primitive.
  }
} else {
  GFq <- seq_len(m)-1                   # Short cut as we do not know the primitive element.
  split(GFq, seq_along(GFq))            # transpose vector
}

q <- length(GFq)                        # Order of field.
stopifnot(q > 1)                        # Non trivial.
if (q != m^d) warning("Polynomial px is not primitive.") # Warning.
if (q %% 4 != 1) warning("Not a Paley graph")            # Warning.

## ----------------------------------------------------------- (iv)
## calculate all quadratic residuals in (QR).
QR <- list()
for ( i in seq_along(GFq[]) ){
    vtx <- GFq[[i]]
    if (BrHa && q == 81) vtx <- gf_mult(vtx, vtx)  # ^2
    QR[[i]] <- gf_mult(vtx, vtx)                   # ^2
}
QR <- unique(QR)

## ------------------------------------------------------------- (v)
## Calculate edges for all vertex pairs in upper half (el)
## Determmine difference k, i.
## If difference is a nonzero quadratic residual then create edge
idx <- 0
df <- data.frame(matrix(ncol = 3, nrow = 0, dimnames=list(NULL, c("to", "from", "diff"))))
for (i in seq_len(q-1)){
  for (k in seq(i+1, q) ){
    idx      <- idx + 1
    diff     <- GFq[[k]] - GFq[[i]]                 # difference between point k, i.
    diff <- (diff + m) %% m                         # all positive.
    stopifnot(!all(diff == 0))                      # nonzero
    df[idx,] <- c(i, k, list(list(diff)) )
  }
}

qresi <- match(df[,3], QR)                          # index in residuals QR.
df[, "qr-idx"] <- qresi
el <- df[!is.na(df[,"qr-idx"]),]                    # Remove non residuals.

## plot and statistics
library(igraph)
g <- graph_from_data_frame(el, directed = FALSE)
k   <- mean(min(degree(g)), max(degree(g)))
g$layout  <- layout_in_circle(g)
g$main <- paste0(ifelse(!exists("BrHa") || BrHa, "Brouwer, Haemers", "Payley"), " graph (v = ", q, ", k = ", k, ")")
V(g)$size <- 5
V(g)$label <- NA
c(q * (q-1) / 4, ecount(g))
plot(g)

## ----------------------------------------------------------- (vi)
ev  <- round(eigen(g[])$values, digits = 6)
nev <- length(unique(ev))
t   <- (q-1) / 4
{
message(sprintf(g$main))
message(sprintf("Vertices               : %i"      , vcount(g)))
message(sprintf("Degree (k)             : %i, %i"  , 2*t/(BrHa+1), k ))
message(sprintf("Eigenvalues (r, s, #)  : %.1f, %.1f, %i", max(ev), min(ev), nev) )
message(sprintf("Diameter               : %i"      , diameter(g)) )
message(sprintf("Radius                 : %i"      , radius(g)) )
message(sprintf("Girth (shortes cycle)  : %i"      , girth(g)$girth) )
message(sprintf("Independent set (α)    : %i"      , length(largest_ivs(g)[[1]])) )
message(sprintf("Edge connectivity      : %.1f"    , edge_connectivity(g)) )
message(sprintf("Edges                  : %.1f, %i", q * t/ (BrHa+1), ecount(g) ))
message(sprintf("Automorphisms          : %.1f, %s", d*q*(q-1) / 2, automorphisms(g)$group_size) )
}
max(greedy_vertex_coloring(g))
stopifnot(min(degree(g)) == max(degree(g)))
stopifnot(max(ev) == max(degree(g)))
stopifnot(nev == 3)

## -------------------------------------------------------------------------------------
## https://www.win.tue.nl/~aeb/graphs/srg/srgtab.html
## SRG parameters.
##  v   - number of vertices
##  k   - valency, degree
##  λ   - number of common neighbours of two adjacent vertices (lambda)
##  μ   - number of common neighbours of two nonadjacent vertices (mu)
##  r^f - positive eigenvalue with multiplicity
##  s^f - negative eigenvalue with multiplicity
##
## Paley graph
## https://www.win.tue.nl/~aeb/graphs/Paley.html
## For q = 4t + 1, the parameters are v = 4t + 1, k = 2t, λ = t − 1, μ = t.
## q 5 9 13 17 25 29 37 41 49 53 61 73 81 89 97 101 109 113 121 125 137 149 157 169 173 181 193 197
## α 2 3  3  3  5  4  4  5  7  5  5  5  9  5  6   5   6   7  11   7   7   7   7  13   8   7   7   8
## χ 3 3  5  6  5  8 10  9  7 11 13 15  9 18 17  21  19  17  11  18  20  22  23  13  22  26  28  25
##
## https://www.win.tue.nl/~aeb/graphs/Brouwer-Haemers.html
## Brouwer Haemers graph, v = 81, k = 20,  λ = 1, μ = 6, α = 15, χ = 6 of 7(E van Dam).
## Radius = 2, Diameter = 2, Girth = 3.
## Automorphisms    : 233,280 (2^6 * 3^6 * 5)
##
## ------------------------------------------------------------- (vii)
## Manual construction of the Paley graph for primes, q = 13.
q <- 13
g13 <- graph_from_literal( 0-1-2-3-4-5-6-7-8-9-10-11-12-0    # (1)
                       , 0-3-6-9-12-2-5-8-11-1-4-7-10-0    # (3)
                       , 0-4-8-12-3-7-11-2-6-10-1-5-9-0    # (4)
     )
dev.new(); plot(g13, layout=layout_in_circle)

# Or more general for q is prime.
q <- 13; QRg <- c(1,3,4)
q <- 17; QRg <- c(1,2,3,5)
q <- 5;  QRg <- c(1)
BrHA <- FALSE
g2 <- graph_from_data_frame(data.frame(from=numeric(), to=numeric()),
                           directed = FALSE,
                           vertices = seq(q) - 1)
for (h in QRg) {
     g2 <- g2 + path(as.character(seq(from=0, to=q*h, by=h) %% q) )
}
dev.new(); plot(g2, layout=layout_in_circle)

## ----------------------------------------------------------------
## compute generators of quadratic residues mod q, q is prime (QR).
q <- 29
q1 <- (seq(q/2) * seq(q/2)) %% q
q2 <- c(q1, -q1 + q) 
QRg <- sort(unique(q2[q2 < q/2]))

[KeesP]

Construction by Blokhuis, Brouwer, Haemers

A more basis approach to create the Brouwer-Haemers graph can be found here:

Blokhuis, A., Brouwer, A. E., & Haemers, W. H. (2012). The graph with spectrum 14^1 2^40 (−4)^10 (−6)^9.Designs, Codes and Cryptography , 65(1-2), 71-75. The graph with spectrum 141 240 (−4)10 (−6)9 | SpringerLink .

It mainly relies on the Hamming distance of two vectors, that is, the number of positions at which the two vectors differ. I could not reconstruct the procedure in the article. We follow a slightly different method.

(i) The first step is to create all the vectors in F35 (not F36 as in the paper), containing all combinations of 0, 1, 2.

library(igraph)
m  <- 3                # K = GF(m), E = GF(m^d).
d  <- 5                # d is degree of extension.
 
## Generate GF(m^d)
F35 <- seq_len(m) - 1L; dfs <- rep(list(F35), d)
F35 <- rev(expand.grid(dfs))            # all combinations over GFm.
q   <- nrow(F35)                        # 3^5.

(ii) The next step is to find the vertices of the graph, the cosets in F35: all vectors with sum is 0 modulo 3.

## ∀x in F35, orthogonal to all ones, c(1,1,1,1,1).
## t(x).<1> == 0 (mod m).
## Add vertices.
coset   <- subset(F35, (rowSums(F35) %% m) == 0)  
cpoints <- apply(coset, 1, paste, collapse="")
n <- nrow(coset)

(iii) Next we construct the edges. Two cosets are adjacent when they differ by a weight-3 vector.

## ------------------------------------------
## Calculate edges for all vertex pairs (el).
## Detemine difference k, i.
## Add edge when they differ by a weight-3 vector.
system.time({
idx <- expand.grid(row = seq_len(n-1), col = seq_len(n)) # all pairs.
idx <- subset(idx, col > row)           # indices upper half matrix.
cd  <- (coset[idx[,1], ] - coset[idx[,2], ] + m) %% m # coset difference.
fi  <- rowSums(cd[,] != 0) == 3          # select weight 3 vectors.
el  <- idx[fi,]                          # edge list Brouwer-Haemers graph.
})

(iv) Create the graph.

g          <- graph_from_data_frame(el, directed = FALSE, vertices=seq(n))
V(g)$name  <- cpoints

Ycoset      <- coset
Yg          <- g
Yg

## IGRAPH 343566e UN-- 81 810 -- 
## + attr: name (v/c)
## + edges from 343566e (vertex names):
##  [1] 00021--00102 00000--00111 00012--00120 00012--00201 00120--00201 00021--00210 00102--00210
##  [8] 00000--00222 00111--00222 00021--01002 00201--01002 00222--01002 00000--01011 00201--01011
## [15] 00210--01011 00012--01020 00210--01020 00222--01020 00000--01101 00021--01101 00120--01101
## [22] 01020--01101 00000--01110 00012--01110 00102--01110 01002--01110 00012--01122 00021--01122
## [29] 00111--01122 01011--01122 00102--01200 00120--01200 00222--01200 01011--01200 01122--01200
## [36] 00102--01212 00111--01212 00201--01212 01020--01212 01101--01212 00111--01221 00120--01221
## [43] 00210--01221 01002--01221 01110--01221 00012--02001 00102--02001 00111--02001 01020--02001
## [50] 01200--02001 01221--02001 00021--02010 00111--02010 00120--02010 01002--02010 01200--02010
## + ... omitted several edges

(v) That is all to it.
Finally, we add some statistics and visualisation .

## ---------------------------------------------------------------------
## (Sub) graph visualizations and statistics
n          <- vcount(g)
k          <- mean(min(degree(g)), max(degree(g)))
g$main     <- paste0("Brouwer-Haemers", " graph (v = ", n, ", k = ", k, ")")
g$sub      <- "Eigenvectors: 1^20, 2^60, (-7)^20"
g$layout   <- layout_in_circle(g)
V(g)$size  <- 5
V(g)$label <- NA
plot(g)

## -----------------------------------------------------------------------
ev  <- round(eigen(g[])$values, digits = 6)
nev <- length(unique(ev))
t   <- (n-1) / 4
S6  <- factorial(6)

{
message(sprintf(g$main))
message(sprintf("Vertices               : %i"      , vcount(g)))
message(sprintf("Degree (k)             : %i", k ))
message(sprintf("Eigenvalues (r, s, #)  : %.1f, %.1f, %i", max(ev), min(ev), nev) )
message(sprintf("Diameter               : %i"      , diameter(g)) )
message(sprintf("Radius                 : %i"      , radius(g)) )
message(sprintf("Girth (shortes cycle)  : %i"      , girth(g)$girth) )
message(sprintf("Independent set (α)    : %i"      , length(largest_ivs(g)[[1]])) )
message(sprintf("Edge connectivity      : %.1f"    , edge_connectivity(g)) )
message(sprintf("Edges                  : %i"      , ecount(g) ))
message(sprintf("Automorphisms          : %s"      , automorphisms(g)$group_size) )
}
rev(table(Eigenvalues=ev))
max(greedy_vertex_coloring(g))
if (min(degree(g)) != max(degree(g))) warning("Not regular, Min, Max degree = ", min(degree(g)), max(degree(g)))
if (nev != 3)                         warning("Not strongly regular, nev = ", nev)

## try some communities
cl <- cluster_louvain(g)
cl <- cluster_leading_eigen(g)
dev.new(); plot(cl, g)

[KeesP]

Subgraphs of the Brouwer-Haemers graphs (Yg)

(i) Generalized Quadrangle GQ(2,1), Rook’s graph.

This graph is the line graph of the complete bipartite graph K3,3.

Construction:

K33  <- make_full_bipartite_graph(3, 3)      # Complete bipartite.
GQ21 <- make_line_graph(K33)                 # Generalized quadrangle, line graph off K33.
vcount(GQ21)
##[1] 9
automorphisms(GQ21)$group_size
##[1] "72"
ev    <- round(eigen(GQ21[])$values, digits = 6)
nev   <- length(unique(ev)); rev(table(Eigenvalues=ev))
## Eigenvalues
## 4  1 -2 
## 1  4  4 
iso21 <- subgraph_isomorphisms(GQ21, Yg, induced=TRUE, time.limit=5 * 60)                        ## Find subgraphs in Yg.
iso21 |> lapply(as_ids) |> lapply(sort) |> sapply(paste, collapse = "-") |> unique() |> length() ## 405 copies.

Visualization:

g              <- GQ21
g$main         <- "Blokhuis-Brouwer-Haemers subgraph\nRook's graph"
g$sub          <- "v = 9, k = 4\nEigenvalues: 4^1  1^4 (-2)^4"
V(g)$name      <- as_ids(iso21[[1]])
E(g)$width     <- 2
E(g)$curved    <- rep(0, ecount(g))
E(g)$curved[c(2, 8, 10, 12, 15, 17)] <- 0.2 
V(g)$label.cex <- 0.8
g$layout       <- layout_on_grid(g)
dev.new(); plot(g)

[KeesP]

(ii) Regular graph Δ of Y, v=60, k=14, ev = 1^14, 2^40, (-4)^10, (-6)^9.

Second subconstituent Δ.

Construction:

## Induced on the set of vertices at distance 2 from a fixed vertex b of Y (81x).
vsg <- V(Yg)[distances(Yg, v = c(1), mode="all") == 2]
Ddg <- subgraph(Yg, v=vsg)                   # Subgraph, vertices v, d(v,b) ≤ 2, for fix b.
setdiff(V(Yg)$name, V(Ddg)$name)
##  [1] "00000" "00111" "00222" "01011" "01101" "01110" "02022" "02202" "02220" "10011"
## [11] "10101" "10110" "11001" "11010" "11100" "20022" "20202" "20220" "22002" "22020"
## [21] "22200"

## Or alternatively, all nonzero Y cosets without: coordinates all 1s, or all 2s.
same <- apply(Ycoset, 1, function(x) length(unique(x[x!=0])) <= 1) == 1
Dsg    <- Yg - V(Yg)[same]                   # Delete vertices.
## sgi    <- subgraph_isomorphisms(Dsg, Yg, induced=TRUE, time.limit=5*60) #automorphism x 3^4
vcount(Dsg)
## [1] 60
ev     <- round(eigen(Dsg[])$values, digits = 6)
nev    <- length(unique(ev)); rev(table(Eigenvalues=ev))
## Eigenvalues
## 14  2 -4 -6 
##  1 40 10  9 

Visualization:

g      <- Dsg; n <- vcount(g)                # Visualisation
g$main <- "Blokhuis-Brouwer-Haemers subgraph (Δ)"
g$sub  <- "Subgraph without vertices with all 0s, 1s, 2s
           v = 60, k = 14
           Eigenvalues: 14^1, 2^40, (-4)^10, (-6)^9"
V(g)$size  <- 5
V(g)$label <- NA
g$layout <- layout_in_circle(g)
dev.new(); plot(g)

[KeesP]

(iii) Root lattices A5 + A5, v = 30, eigenvalues: 8^1, 2^15, (-2)^9 (-4)^5.

Small regular graphs with four eigenvalues, Appendix A No. 130.

Construction:

same <- apply(Ycoset, 1, function(x) length(unique(x[x!=0])) <= 1) == 1 # from Yg
single_zero   <- rowSums(Ycoset[,] == 0) == 1
A5ga          <- subgraph(Yg, V(Yg)[!same &  single_zero]) # First  summand
A5gb          <- subgraph(Yg, V(Yg)[!same & !single_zero]) # Second summand
subgraph_isomorphic(A5gb, Dsg, induced=TRUE)
## [1] TRUE
subgraph_isomorphic(A5ga, Dsg, induced=TRUE)
## [1] TRUE
vcount(A5ga)
[1] 30
ev     <- round(eigen(A5ga[])$values, digits = 6)
nev    <- length(unique(ev)); rev(table(Eigenvalues=ev))
## Eigenvalues
## 8  2 -2 -4 
## 1 15  9  5 

Visualization:

g <- A5ga
g$main        <- "Blokhuis-Brouwer-Haemers subgraph (Δ = A5 + A5)\nRoot lattice A5"
g$sub         <- "Subgraph with (or without) single zero in GF(3^5)\nv = 30, k = 8
                 Eigenvalues: 8^1, 2^15, (-2)^9, (-4)^5"
V(g)$label.cex <- 0.6
E(g)$width     <- 2
lll <- layout_on_grid(g, width=6)
dev.new(); plot(g, layout=lll)

Note:
I was not able to construct a nice symmetric layout of A5 as in Root system (Wikipedia) .

.