Example: Rounding the Delsarte LP bound

For this example, we slightly modify the code of the Delsarte bound to define an exact problem.

using ClusteredLowRankSolver, Nemo

function delsarte_exact(n, d, costheta; FF=QQ, g=1, eps=1e-40)
    constraints = []

    P, x = polynomial_ring(FF, :x)

    gbasis = basis_gegenbauer(2d, n, x)
    sosbasis = basis_chebyshev(2d, x)

    samples = sample_points_chebyshev(2d)
    # round the samples to QQ:
    samples = [round(BigInt, x * 10^4)//10^4 for x in samples]

    c = Dict()
    for k = 0:2d
        c[k] = [gbasis[k+1];;]
    c[:A] = LowRankMatPol([1], [sosbasis[1:d+1]])
    c[:B] = LowRankMatPol([(x+1)*(costheta-x)], [sosbasis[1:d]])
    push!(constraints, Constraint(-1, c, Dict(), samples))

    objective = Objective(1, Dict(k => [1;;] for k=0:2d), Dict())

    problem = Problem(Minimize(objective), constraints)

    problem_bigfloat = map(x->generic_embedding(x, g), problem)
    status, primalsol, dualsol, time, errorcode = solvesdp(problem_bigfloat, duality_gap_threshold=eps)

    return objvalue(problem, dualsol), problem, primalsol, dualsol

Here we made a few modifications.

  • We defined the polynomial ring over the field FF, so that the problem becomes exact.
  • We also made the samples exact rational numbers. This is in this case not strictly necessary if we use coefficient matching for the rounding procedure, but it is necessary when evaluating polynomials on these samples.
  • We added the line
        problem_bigfloat = map(x->generic_embedding(x, g), problem)
    This embeds the field in $\mathbb{R}$ using the floating point approximation g of the generator of the field we use. The exact problem cannot be solved with the solver if we use a number field with a generator, since it is unclear which generator of the field we meant when building the problem.
  • We now return the problem and the primal and dual solution too.

One example where the delsarte bound is sharp, is when considering a code with $\cos\theta = 1/(\sqrt{5} - 1)$ in $\mathbb{R}^4$. The optimal spherical code then has $120$ points. To round the solution, we first define the field using the minimal polynomial $x^2 - 5 = 0$.

R, x = polynomial_ring(QQ, :x)
N, z = number_field(x^2 - 5, :z)
gapprox = sqrt(big(5))

Then we can call the rounding procedure:

n = 4
d = 9
costheta = 1/(z-1)
# find an approximate solution
obj, problem, primalsol, dualsol = delsarte_exact(n, d, costheta; FF=N, g = gapprox)
# round the approximate solution to an exact solution
success, exactdualsol = exact_solution(problem, primalsol, dualsol; FF=N, g=gapprox)
objexact = objvalue(problem, exactdualsol)
(success, objexact)
(true, 120)

When the problem is not defined over the same field as the solution, we can find the field using

N2, gapprox2 = find_field(primalsol, dualsol)
defining_polynomial(N2), gapprox2
(x^2 + x - 1, -1.618033988749894848204586834365638117720309179805762862135448622705260462818892)

which returns the field with defining polynomial $x^2 + x - 1 = 0$ and approximate generator $-1.618033... \approx (-1 -\sqrt{5})/2$.

Using coefficient matching

As mentioned in the section about Rounding, using coefficient matching will often result in a solution of smaller bit size. To do that in this example, one can use

obj, problem, primalsol, dualsol = delsarte_exact(n, d, costheta; FF=N, g = gapprox)
R, x = polynomial_ring(N, :x)
mon_basis = [x^k for k=0:2d]
success, exactdualsol_smaller = exact_solution(problem, primalsol, dualsol, FF=N, g=gapprox, monomial_bases = [mon_basis])
 iter  time(s)           μ       P-obj       D-obj        gap    P-error    p-error    d-error        α_p        α_d       beta
    1      0.1   1.000e+20   1.000e+00   1.900e+11   1.00e+00   1.00e+10   0.00e+00   2.18e+11   3.69e-01   5.95e-01   3.00e-01
    2      0.1   6.494e+19   1.223e+10   1.739e+11   8.69e-01   6.31e+09   0.00e+00   8.84e+10   7.31e-01   6.03e-01   3.00e-01
    3      0.1   2.817e+19   3.102e+10   2.208e+11   7.54e-01   1.70e+09   0.00e+00   3.51e+10   6.85e-01   7.10e-01   3.00e-01
    4      0.1   1.230e+19   3.546e+10   3.600e+11   8.21e-01   5.34e+08   0.00e+00   1.02e+10   5.57e-01   1.00e+00   3.00e-01
    5      0.1   8.216e+18   2.178e+10   8.065e+11   9.47e-01   2.37e+08   0.00e+00   6.38e-65   7.69e-01   1.00e+00   3.00e-01
    6      0.1   3.035e+18   5.560e+09   1.290e+12   9.91e-01   5.47e+07   0.00e+00   5.00e-64   8.01e-01   1.00e+00   3.00e-01
    7      0.1   9.665e+17   1.150e+09   2.064e+12   9.99e-01   1.09e+07   0.00e+00   7.96e-64   8.65e-01   1.00e+00   3.00e-01
    8      0.1   2.092e+17   1.573e+08   3.302e+12   1.00e+00   1.47e+06   0.00e+00   3.18e-63   8.98e-01   1.00e+00   3.00e-01
    9      0.1   3.428e+16   1.603e+07   5.284e+12   1.00e+00   1.51e+05   0.00e+00   1.40e-63   8.88e-01   1.00e+00   3.00e-01
   10      0.2   6.127e+15   1.797e+06   8.453e+12   1.00e+00   1.68e+04   0.00e+00   2.62e-63   8.99e-01   1.00e+00   3.00e-01
   11      0.2   9.935e+14   1.816e+05   1.352e+13   1.00e+00   1.71e+03   0.00e+00   1.15e-63   8.93e-01   1.00e+00   3.00e-01
   12      0.2   1.699e+14   1.946e+04   2.163e+13   1.00e+00   1.82e+02   0.00e+00   2.88e-63   9.00e-01   1.00e+00   3.00e-01
   13      0.2   2.794e+13   2.009e+03   3.442e+13   1.00e+00   1.82e+01   0.00e+00   1.79e-62   8.98e-01   1.00e+00   3.00e-01
   14      0.2   5.597e+12   2.662e+02   5.231e+13   1.00e+00   1.86e+00   0.00e+00   2.02e-62   8.79e-01   1.00e+00   3.00e-01
   15      0.2   2.030e+12   9.171e+01   5.562e+13   1.00e+00   2.25e-01   0.00e+00   5.15e-63   7.97e-01   1.00e+00   3.00e-01
   16      0.2   7.056e+11   7.350e+01   2.417e+13   1.00e+00   4.58e-02   0.00e+00   2.40e-62   8.24e-01   1.00e+00   3.00e-01
   17      0.2   2.136e+11   7.073e+01   7.703e+12   1.00e+00   8.06e-03   0.00e+00   1.40e-63   1.00e+00   1.00e+00   3.00e-01
   18      0.2   6.305e+10   6.979e+01   2.396e+12   1.00e+00   1.11e-75   0.00e+00   1.74e-62   1.00e+00   1.00e+00   3.00e-01
   19      0.3   1.891e+10   6.985e+01   7.188e+11   1.00e+00   5.53e-76   0.00e+00   1.01e-61   9.94e-01   9.94e-01   1.00e-01
   20      0.3   1.996e+09   6.986e+01   7.583e+10   1.00e+00   1.11e-75   0.00e+00   3.97e-63   1.00e+00   1.00e+00   1.00e-01
   21      0.3   2.003e+08   6.986e+01   7.613e+09   1.00e+00   1.11e-75   0.00e+00   3.60e-64   1.00e+00   1.00e+00   1.00e-01
   22      0.3   2.005e+07   6.987e+01   7.619e+08   1.00e+00   1.11e-75   0.00e+00   1.70e-65   1.00e+00   1.00e+00   1.00e-01
   23      0.3   2.005e+06   6.987e+01   7.619e+07   1.00e+00   5.53e-76   0.00e+00   1.35e-66   1.00e+00   1.00e+00   1.00e-01
   24      0.3   2.005e+05   6.988e+01   7.620e+06   1.00e+00   1.11e-75   0.00e+00   2.72e-67   1.00e+00   1.00e+00   1.00e-01
   25      0.3   2.006e+04   6.988e+01   7.622e+05   1.00e+00   1.11e-75   0.00e+00   8.00e-69   1.00e+00   1.00e+00   1.00e-01
   26      0.3   2.007e+03   6.989e+01   7.635e+04   9.98e-01   1.11e-75   0.00e+00   1.69e-69   9.99e-01   9.99e-01   1.00e-01
   27      0.3   2.026e+02   6.998e+01   7.768e+03   9.82e-01   1.11e-75   0.00e+00   4.16e-70   9.90e-01   9.90e-01   1.00e-01
   28      0.3   2.205e+01   7.086e+01   9.087e+02   8.55e-01   1.11e-75   0.00e+00   4.92e-71   9.26e-01   9.26e-01   1.00e-01
   29      0.4   3.666e+00   7.788e+01   2.172e+02   4.72e-01   1.11e-75   0.00e+00   2.12e-71   8.10e-01   8.10e-01   1.00e-01
   30      0.4   9.925e-01   1.015e+02   1.392e+02   1.57e-01   1.11e-75   0.00e+00   2.28e-71   6.72e-01   6.72e-01   1.00e-01
   31      0.4   3.920e-01   1.120e+02   1.269e+02   6.23e-02   5.53e-76   0.00e+00   6.82e-71   8.04e-01   8.04e-01   1.00e-01
   32      0.4   1.082e-01   1.179e+02   1.220e+02   1.71e-02   5.53e-76   0.00e+00   5.75e-71   8.72e-01   8.72e-01   1.00e-01
   33      0.4   2.330e-02   1.195e+02   1.204e+02   3.69e-03   1.11e-75   0.00e+00   4.80e-71   9.67e-01   9.67e-01   1.00e-01
   34      0.4   3.027e-03   1.199e+02   1.201e+02   4.79e-04   1.11e-75   0.00e+00   5.06e-71   9.83e-01   9.83e-01   1.00e-01
   35      0.4   3.477e-04   1.200e+02   1.200e+02   5.51e-05   1.11e-75   0.00e+00   3.86e-71   9.94e-01   9.94e-01   1.00e-01
   36      0.4   3.680e-05   1.200e+02   1.200e+02   5.83e-06   3.32e-75   0.00e+00   9.91e-72   9.99e-01   9.99e-01   1.00e-01
   37      0.4   3.724e-06   1.200e+02   1.200e+02   5.90e-07   3.32e-75   0.00e+00   3.08e-71   1.00e+00   1.00e+00   1.00e-01
   38      0.5   3.730e-07   1.200e+02   1.200e+02   5.91e-08   1.11e-75   0.00e+00   1.42e-71   1.00e+00   1.00e+00   1.00e-01
   39      0.5   3.731e-08   1.200e+02   1.200e+02   5.91e-09   1.11e-75   0.00e+00   2.61e-71   1.00e+00   1.00e+00   1.00e-01
   40      0.5   3.732e-09   1.200e+02   1.200e+02   5.91e-10   1.11e-75   0.00e+00   4.04e-71   1.00e+00   1.00e+00   1.00e-01
   41      0.5   3.732e-10   1.200e+02   1.200e+02   5.91e-11   2.21e-75   0.00e+00   7.64e-71   1.00e+00   1.00e+00   1.00e-01
   42      0.5   3.733e-11   1.200e+02   1.200e+02   5.91e-12   1.11e-75   0.00e+00   6.53e-72   1.00e+00   1.00e+00   1.00e-01
   43      0.5   3.733e-12   1.200e+02   1.200e+02   5.91e-13   1.11e-75   0.00e+00   3.66e-71   1.00e+00   1.00e+00   1.00e-01
   44      0.5   3.733e-13   1.200e+02   1.200e+02   5.91e-14   1.11e-75   0.00e+00   5.06e-71   1.00e+00   1.00e+00   1.00e-01
   45      0.5   3.734e-14   1.200e+02   1.200e+02   5.91e-15   5.53e-76   0.00e+00   1.73e-70   1.00e+00   1.00e+00   1.00e-01
   46      0.5   3.734e-15   1.200e+02   1.200e+02   5.91e-16   1.11e-75   0.00e+00   8.91e-71   1.00e+00   1.00e+00   1.00e-01
   47      0.5   3.734e-16   1.200e+02   1.200e+02   5.91e-17   1.11e-75   0.00e+00   1.03e-70   1.00e+00   1.00e+00   1.00e-01
   48      0.6   3.735e-17   1.200e+02   1.200e+02   5.91e-18   5.53e-76   0.00e+00   2.09e-70   1.00e+00   1.00e+00   1.00e-01
   49      0.6   3.735e-18   1.200e+02   1.200e+02   5.91e-19   2.21e-75   0.00e+00   8.53e-70   1.00e+00   1.00e+00   1.00e-01
   50      0.6   3.736e-19   1.200e+02   1.200e+02   5.91e-20   3.32e-75   0.00e+00   2.10e-69   1.00e+00   1.00e+00   1.00e-01
   51      0.6   3.736e-20   1.200e+02   1.200e+02   5.92e-21   1.11e-75   0.00e+00   1.44e-69   1.00e+00   1.00e+00   1.00e-01
   52      0.6   3.736e-21   1.200e+02   1.200e+02   5.92e-22   1.11e-75   0.00e+00   6.31e-69   1.00e+00   1.00e+00   1.00e-01
   53      0.6   3.737e-22   1.200e+02   1.200e+02   5.92e-23   1.11e-75   0.00e+00   7.83e-69   1.00e+00   1.00e+00   1.00e-01
   54      0.6   3.737e-23   1.200e+02   1.200e+02   5.92e-24   1.11e-75   0.00e+00   6.78e-69   1.00e+00   1.00e+00   1.00e-01
   55      0.6   3.737e-24   1.200e+02   1.200e+02   5.92e-25   5.53e-76   0.00e+00   8.59e-68   1.00e+00   1.00e+00   1.00e-01
   56      0.6   3.738e-25   1.200e+02   1.200e+02   5.92e-26   1.11e-75   0.00e+00   1.01e-67   1.00e+00   1.00e+00   1.00e-01
   57      0.7   3.738e-26   1.200e+02   1.200e+02   5.92e-27   2.21e-75   0.00e+00   7.42e-68   1.00e+00   1.00e+00   1.00e-01
   58      0.7   3.739e-27   1.200e+02   1.200e+02   5.92e-28   2.21e-75   0.00e+00   1.44e-66   1.00e+00   1.00e+00   1.00e-01
   59      0.7   3.739e-28   1.200e+02   1.200e+02   5.92e-29   1.11e-75   0.00e+00   1.19e-66   1.00e+00   1.00e+00   1.00e-01
   60      0.7   3.739e-29   1.200e+02   1.200e+02   5.92e-30   5.53e-76   0.00e+00   5.83e-66   1.00e+00   1.00e+00   1.00e-01
   61      0.7   3.740e-30   1.200e+02   1.200e+02   5.92e-31   1.11e-75   0.00e+00   6.90e-66   1.00e+00   1.00e+00   1.00e-01
   62      0.7   3.740e-31   1.200e+02   1.200e+02   5.92e-32   1.11e-75   0.00e+00   9.38e-66   1.00e+00   1.00e+00   1.00e-01
   63      0.7   3.740e-32   1.200e+02   1.200e+02   5.92e-33   1.11e-75   0.00e+00   2.37e-66   1.00e+00   1.00e+00   1.00e-01
   64      0.7   3.741e-33   1.200e+02   1.200e+02   5.92e-34   1.11e-75   0.00e+00   1.26e-64   1.00e+00   1.00e+00   1.00e-01
   65      0.7   3.741e-34   1.200e+02   1.200e+02   5.92e-35   2.21e-75   0.00e+00   7.82e-65   1.00e+00   1.00e+00   1.00e-01
   66      0.7   3.742e-35   1.200e+02   1.200e+02   5.92e-36   1.11e-75   0.00e+00   6.95e-65   1.00e+00   1.00e+00   1.00e-01
   67      0.8   3.742e-36   1.200e+02   1.200e+02   5.92e-37   1.11e-75   0.00e+00   3.23e-64   1.00e+00   1.00e+00   1.00e-01
   68      0.8   3.742e-37   1.200e+02   1.200e+02   5.93e-38   2.21e-75   0.00e+00   1.24e-63   1.00e+00   1.00e+00   1.00e-01
   69      0.8   3.743e-38   1.200e+02   1.200e+02   5.93e-39   1.11e-75   0.00e+00   3.03e-63   1.00e+00   1.00e+00   1.00e-01
   70      0.8   3.743e-39   1.200e+02   1.200e+02   5.93e-40   1.11e-75   0.00e+00   6.50e-63   1.00e+00   1.00e+00   1.00e-01
Optimal solution found
  0.796961 seconds (1.53 M allocations: 118.700 MiB, 15.43% gc time)
 iter  time(s)           μ       P-obj       D-obj        gap    P-error    p-error    d-error        α_p        α_d       beta

Primal objective:119.99999999999999999999999999999999999999176437133780968081009238851823100077
Dual objective:120.0000000000000000000000000000000000000059895481179565957744782628958319998673
Duality gap:5.927156991727881235160780990667082957200110525181464317207184020991804734467005e-41
** Starting computation of basis transformations **
  Block 14 of size 1 x 1
  Block 3 of size 1 x 1
  Block 17 of size 1 x 1
  Block 6 of size 1 x 1
  Block 9 of size 1 x 1
  Block 12 of size 1 x 1
    Block 12 has 1 kernel vectors. The maximum numerator and denominator are 1 and 1
    After reduction, the maximum number of the basis transformation matrix is 1
  Block 1 of size 1 x 1
  Block 15 of size 1 x 1
  Block 4 of size 1 x 1
  Block 18 of size 1 x 1
  Block 0 of size 1 x 1
    Block 0 has 1 kernel vectors. The maximum numerator and denominator are 1 and 1
    After reduction, the maximum number of the basis transformation matrix is 1
  Block 7 of size 1 x 1
  Block 10 of size 1 x 1
  Block 13 of size 1 x 1
  Block 2 of size 1 x 1
  Block 16 of size 1 x 1
  Block 5 of size 1 x 1
  Block 8 of size 1 x 1
  Block 11 of size 1 x 1
  Block B of size 9 x 9
    Block B has 6 kernel vectors. The maximum numerator and denominator are 18 and 2
    After reduction, the maximum number of the basis transformation matrix is 10
  Block A of size 10 x 10
    Block A has 8 kernel vectors. The maximum numerator and denominator are 12 and 1
    After reduction, the maximum number of the basis transformation matrix is 1
** Finished computation of basis transformations (0.016241973s) **
** Transforming the problem and the solution ** (0.003062119s)
** Projection the solution into the affine space **
  Reducing the system from 26 columns to 26 columns
  Constructing the linear system... (0.106886787s)
  Computing an approximate solution in the extension field... (0.004195144s)
  Preprocessing to get an integer system... (0.000846855s)
  Finding the pivots of A using RREF mod p... (0.001758689 0.047980796 s)
  Solving the system of size 38 x 40 using the pseudoinverse... 0.002755552s
** Finished projection into affine space (0.167724697s) **
** Checking feasibility **
The slacks are satisfied (checked or ensured by solving the system)
Checking sdp constraints
 done (0.020462578)

This is recommended especially for larger polynomial programs.

Getting an exact feasible solution close to the optimum

To modify the code to find an exact strictly feasible solution, we change the code for the objective to the following:

   if isnothing(obj)
      objective = Objective(1, Dict(k => [1;;] for k=0:2d), Dict())
      objective = Objective(0, Dict(), Dict()) 
      push!(constraints, Constraint(obj-1, Dict(k => [1;;] for k=0:2d), Dict()))

where the function signature is now given by

function delsarte_exact(n, d, costheta; obj=nothing, FF=QQ, g=1, eps=1e-40)

That is, it takes an extra argument obj. When obj is set to a numerical value, this adds the constraint that the objective should be equal to obj. Then the exact solution can be found as follows.

d = 10
# find the objective
obj_initial, problem_initial, _, _ = delsarte_exact(3, d, 1//2)
# find a strictly feasible solution with a slightly larger objective
obj = obj_initial + 1e-6
_, problem, primalsol, dualsol = delsarte_exact(3, d, 1//2, obj=rationalize(obj))
# round the solution
successfull, exactdualsol = exact_solution(problem, primalsol, dualsol)
Float64(objvalue(problem_initial, exactdualsol)), obj_initial
(13.158315347390305, 13.15831434739030571916973180897335312215428948322487489521314832155814709836958)