GitHub

The Theta'-function of Erdos-Renyi graphs

Let $q$ be an odd prime, and let $V = \mathrm{GF}(q)^3$ be a three dimensional vector space over the finite field of order $q$. The set of one dimensional subspaces, i.e. the projective plane, of $V$ is denoted by $\mathrm{PG}(2,q)$. There are $q^2+q+1$ such subspaces, which are the vertices of the Erdos-Renyi graph $\mathrm{ER}(q)$. Two vertices are adjacent if they are distinct and orthogonal, i.e. for two representing vectors $x$ and $y$ we have $x^Ty=0$. We are interested in the size of a maximum stable set of these graphs, specifically upper bounds for this value. Note that these are not the equally named random graphs.

q = 7
PG2q = vcat([[0, 0, 1]],
    [[0, 1, b] for b = 0:q-1],
    [[1, a, b] for a = 0:q-1 for b = 0:q-1])
Adj = [x' * y % q == 0 && x != y for x in PG2q, y in PG2q]
size(Adj)
(57, 57)
spy(Adj)
Example block output

The Theta'-function

The Theta'-function $\vartheta(G)$ of a graph $G=(V,E)$ is such an upper bound, based on semidefinite programming:

\[\vartheta(G)\coloneqq \sup\{\langle X,J\rangle : \langle X,A\rangle = 0, X\succcurlyeq 0, X\geq 0\}.\]

In vectorized standard form this is simply

N = length(PG2q)
C = ones(N^2)
A = vcat(vec(Adj)', vec(Matrix{Float64}(I, N, N))')
b = [0.0, 1.0];

Determining the symmetry reduction

We can now apply the Jordan reduction method to the problem. First, we need to determine an (optimal) admissible subspace.

using SDPSymmetryReduction
P = admissible_subspace(C, A, b; verbose=true)
dim(P)

using Test #src
┌ Info: Starting the reduction. Dimensions:
│   maximal = 1653
└   initial = 2
┌ Info: Minimal admissible subspace converged in 5 iterations at dimension:
└   final = 18

Running admissible_subspace returns a Partition object. dim(P) are the number of orbits (and thus variables), and P.matrix is a matrix with integer values from 1 trough dim(P). Here, P.matrix looks like this (different color shades = different orbits):

Example block output

Now we can block-diagonalize the algebra (numerically)

blkD = blockDiagonalize(P, true);
[ Info: Determining eigen-decomposition over Float64...
[ Info:   0.001199 seconds (67 allocations: 153.250 KiB)
[ Info: Determining the algebra-isomorphism...
[ Info:   0.000061 seconds (120 allocations: 56.234 KiB)
[ Info: Calculating image of the basis of the algebra...
[ Info:   0.000237 seconds (362 allocations: 68.203 KiB)

Building the reduced SDP

Since blkD.blks[i] is the block-diagonalized image of P.matrix .== i, we obtain the new, symmetry reduced SDP by

using SparseArrays
PMat = hcat([sparse(vec(P.matrix .== i)) for i = 1:dim(P)]...)
newA = A * PMat
newB = b
newC = C' * PMat;

Solving the SDP with JuMP and CSDP

using JuMP, CSDP
m = Model(CSDP.Optimizer)

# Initialize variables corresponding parts of the partition P
# >= 0 because the original SDP-matrices are entry-wise nonnegative
x = @variable(m, x[1:dim(P)] >= 0)

@constraint(m, newA * x .== newB)
@objective(m, Max, newC * x)

psdBlocks = sum(blkD.blks[i] .* x[i] for i = 1:dim(P))
for blk in psdBlocks
    if size(blk, 1) > 1
        @constraint(m, blk in PSDCone())
    else
        @constraint(m, blk .>= 0)
    end
end

optimize!(m)
Iter: 21 Ap: 9.57e-01 Pobj: -7.7942186e+00 Ad: 9.56e-01 Dobj: -7.7942189e+00
Success: SDP solved
Primal objective value: -7.7942186e+00
Dual objective value: -7.7942189e+00
Relative primal infeasibility: 7.02e-13
Relative dual infeasibility: 4.86e-10
Real Relative Gap: -1.30e-08
XZ Relative Gap: 9.99e-10
DIMACS error measures: 7.03e-13 0.00e+00 1.80e-08 0.00e+00 -1.30e-08 9.99e-10
CSDP 6.2.0
Iter:  0 Ap: 0.00e+00 Pobj:  3.2490000e+04 Ad: 0.00e+00 Dobj:  0.0000000e+00
Iter:  1 Ap: 8.64e-02 Pobj:  3.4910670e+04 Ad: 3.84e-01 Dobj:  8.5393179e+00
Iter:  2 Ap: 3.98e-01 Pobj:  1.7915117e+04 Ad: 1.00e+00 Dobj:  4.6958195e+01
Iter:  3 Ap: 8.92e-01 Pobj:  1.9764078e+03 Ad: 1.00e+00 Dobj:  4.9990490e+01
Iter:  4 Ap: 9.50e-01 Pobj:  1.0706885e+02 Ad: 1.00e+00 Dobj:  5.0266844e+01
Iter:  5 Ap: 9.21e-01 Pobj:  1.7107795e+01 Ad: 1.00e+00 Dobj:  3.3879740e+01
Iter:  6 Ap: 8.54e-01 Pobj:  1.4377274e+01 Ad: 9.82e-01 Dobj:  2.2336813e+01
Iter:  7 Ap: 7.22e-01 Pobj:  1.4872793e+01 Ad: 1.00e+00 Dobj:  1.7272339e+01
Iter:  8 Ap: 8.41e-01 Pobj:  1.5516454e+01 Ad: 1.00e+00 Dobj:  1.6183691e+01
Iter:  9 Ap: 8.76e-01 Pobj:  1.5676799e+01 Ad: 1.00e+00 Dobj:  1.5879167e+01
Iter: 10 Ap: 7.83e-01 Pobj:  1.5720533e+01 Ad: 1.00e+00 Dobj:  1.5765936e+01
Iter: 11 Ap: 9.11e-01 Pobj:  1.5740504e+01 Ad: 1.00e+00 Dobj:  1.5745659e+01
Iter: 12 Ap: 9.29e-01 Pobj:  1.5743131e+01 Ad: 1.00e+00 Dobj:  1.5743456e+01
Iter: 13 Ap: 9.49e-01 Pobj:  1.5743384e+01 Ad: 1.00e+00 Dobj:  1.5743290e+01
Iter: 14 Ap: 9.42e-01 Pobj:  1.5743401e+01 Ad: 1.00e+00 Dobj:  1.5743400e+01
termination_status(m)
OPTIMAL::TerminationStatusCode = 1
objective_value(m)
15.74340268223059

This page was generated using Literate.jl.