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];;]
end
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
endHere 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
This embeds the field in $\mathbb{R}$ using the floating point approximationproblem_bigfloat = map(x->generic_embedding(x, g), problem)gof 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.0 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 7.77e-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 6.07e-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 6.19e-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 2.13e-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.85e-63 8.88e-01 1.00e+00 3.00e-01
10 0.1 6.127e+15 1.797e+06 8.453e+12 1.00e+00 1.68e+04 0.00e+00 5.46e-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.63e-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 1.10e-62 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 6.52e-63 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 3.65e-62 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 1.96e-63 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 4.05e-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 6.38e-62 1.00e+00 1.00e+00 3.00e-01
19 0.2 1.891e+10 6.985e+01 7.188e+11 1.00e+00 1.11e-75 0.00e+00 3.22e-61 9.94e-01 9.94e-01 1.00e-01
20 0.2 1.996e+09 6.986e+01 7.583e+10 1.00e+00 2.21e-75 0.00e+00 2.83e-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 5.53e-76 0.00e+00 2.09e-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 2.11e-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 1.11e-75 0.00e+00 5.81e-67 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 1.36e-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 4.41e-69 1.00e+00 1.00e+00 1.00e-01
26 0.3 2.008e+03 6.989e+01 7.636e+04 9.98e-01 5.53e-76 0.00e+00 3.36e-70 9.99e-01 9.99e-01 1.00e-01
27 0.3 2.026e+02 6.998e+01 7.769e+03 9.82e-01 1.11e-75 0.00e+00 7.74e-71 9.90e-01 9.90e-01 1.00e-01
28 0.3 2.205e+01 7.086e+01 9.088e+02 8.55e-01 1.11e-75 0.00e+00 1.19e-70 9.26e-01 9.26e-01 1.00e-01
29 0.3 3.667e+00 7.788e+01 2.172e+02 4.72e-01 1.11e-75 0.00e+00 1.92e-70 8.10e-01 8.10e-01 1.00e-01
30 0.3 9.926e-01 1.015e+02 1.392e+02 1.57e-01 5.53e-76 0.00e+00 3.17e-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 1.11e-75 0.00e+00 9.60e-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 1.11e-75 0.00e+00 8.04e-71 8.72e-01 8.72e-01 1.00e-01
33 0.4 2.331e-02 1.195e+02 1.204e+02 3.69e-03 1.11e-75 0.00e+00 8.21e-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 3.80e-71 9.83e-01 9.83e-01 1.00e-01
35 0.4 3.478e-04 1.200e+02 1.200e+02 5.51e-05 1.11e-75 0.00e+00 2.01e-71 9.94e-01 9.94e-01 1.00e-01
36 0.4 3.681e-05 1.200e+02 1.200e+02 5.83e-06 1.11e-75 0.00e+00 1.58e-71 9.99e-01 9.99e-01 1.00e-01
37 0.4 3.725e-06 1.200e+02 1.200e+02 5.90e-07 2.21e-75 0.00e+00 2.31e-71 1.00e+00 1.00e+00 1.00e-01
38 0.4 3.731e-07 1.200e+02 1.200e+02 5.91e-08 5.53e-76 0.00e+00 3.43e-71 1.00e+00 1.00e+00 1.00e-01
39 0.4 3.732e-08 1.200e+02 1.200e+02 5.91e-09 1.11e-75 0.00e+00 4.59e-71 1.00e+00 1.00e+00 1.00e-01
40 0.4 3.733e-09 1.200e+02 1.200e+02 5.91e-10 5.53e-76 0.00e+00 1.04e-71 1.00e+00 1.00e+00 1.00e-01
41 0.5 3.733e-10 1.200e+02 1.200e+02 5.91e-11 1.11e-75 0.00e+00 1.67e-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 2.59e-71 1.00e+00 1.00e+00 1.00e-01
43 0.5 3.734e-12 1.200e+02 1.200e+02 5.91e-13 1.11e-75 0.00e+00 1.05e-70 1.00e+00 1.00e+00 1.00e-01
44 0.5 3.734e-13 1.200e+02 1.200e+02 5.91e-14 2.21e-75 0.00e+00 7.67e-71 1.00e+00 1.00e+00 1.00e-01
45 0.5 3.735e-14 1.200e+02 1.200e+02 5.91e-15 2.21e-75 0.00e+00 6.24e-71 1.00e+00 1.00e+00 1.00e-01
46 0.5 3.735e-15 1.200e+02 1.200e+02 5.91e-16 1.11e-75 0.00e+00 3.05e-70 1.00e+00 1.00e+00 1.00e-01
47 0.5 3.735e-16 1.200e+02 1.200e+02 5.91e-17 1.11e-75 0.00e+00 4.78e-70 1.00e+00 1.00e+00 1.00e-01
48 0.5 3.736e-17 1.200e+02 1.200e+02 5.91e-18 1.11e-75 0.00e+00 7.27e-70 1.00e+00 1.00e+00 1.00e-01
49 0.5 3.736e-18 1.200e+02 1.200e+02 5.92e-19 1.11e-75 0.00e+00 2.53e-69 1.00e+00 1.00e+00 1.00e-01
50 0.5 3.736e-19 1.200e+02 1.200e+02 5.92e-20 1.11e-75 0.00e+00 6.88e-70 1.00e+00 1.00e+00 1.00e-01
51 0.6 3.737e-20 1.200e+02 1.200e+02 5.92e-21 1.11e-75 0.00e+00 1.25e-68 1.00e+00 1.00e+00 1.00e-01
52 0.6 3.737e-21 1.200e+02 1.200e+02 5.92e-22 1.11e-75 0.00e+00 1.14e-68 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 4.53e-68 1.00e+00 1.00e+00 1.00e-01
54 0.6 3.738e-23 1.200e+02 1.200e+02 5.92e-24 1.11e-75 0.00e+00 1.22e-68 1.00e+00 1.00e+00 1.00e-01
55 0.6 3.738e-24 1.200e+02 1.200e+02 5.92e-25 1.11e-75 0.00e+00 1.80e-67 1.00e+00 1.00e+00 1.00e-01
56 0.6 3.739e-25 1.200e+02 1.200e+02 5.92e-26 1.11e-75 0.00e+00 1.39e-67 1.00e+00 1.00e+00 1.00e-01
57 0.6 3.739e-26 1.200e+02 1.200e+02 5.92e-27 1.11e-75 0.00e+00 1.53e-67 1.00e+00 1.00e+00 1.00e-01
58 0.6 3.739e-27 1.200e+02 1.200e+02 5.92e-28 2.21e-75 0.00e+00 1.03e-67 1.00e+00 1.00e+00 1.00e-01
59 0.6 3.740e-28 1.200e+02 1.200e+02 5.92e-29 3.32e-75 0.00e+00 6.27e-66 1.00e+00 1.00e+00 1.00e-01
60 0.6 3.740e-29 1.200e+02 1.200e+02 5.92e-30 1.11e-75 0.00e+00 1.02e-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 9.50e-66 1.00e+00 1.00e+00 1.00e-01
62 0.7 3.741e-31 1.200e+02 1.200e+02 5.92e-32 1.11e-75 0.00e+00 4.01e-65 1.00e+00 1.00e+00 1.00e-01
63 0.7 3.741e-32 1.200e+02 1.200e+02 5.92e-33 1.11e-75 0.00e+00 2.81e-65 1.00e+00 1.00e+00 1.00e-01
64 0.7 3.742e-33 1.200e+02 1.200e+02 5.92e-34 2.21e-75 0.00e+00 1.01e-65 1.00e+00 1.00e+00 1.00e-01
65 0.7 3.742e-34 1.200e+02 1.200e+02 5.92e-35 1.11e-75 0.00e+00 1.81e-64 1.00e+00 1.00e+00 1.00e-01
66 0.7 3.742e-35 1.200e+02 1.200e+02 5.93e-36 1.11e-75 0.00e+00 2.61e-64 1.00e+00 1.00e+00 1.00e-01
67 0.7 3.743e-36 1.200e+02 1.200e+02 5.93e-37 1.11e-75 0.00e+00 7.49e-65 1.00e+00 1.00e+00 1.00e-01
68 0.7 3.743e-37 1.200e+02 1.200e+02 5.93e-38 1.11e-75 0.00e+00 1.15e-64 1.00e+00 1.00e+00 1.00e-01
69 0.7 3.743e-38 1.200e+02 1.200e+02 5.93e-39 1.11e-75 0.00e+00 1.32e-63 1.00e+00 1.00e+00 1.00e-01
70 0.7 3.744e-39 1.200e+02 1.200e+02 5.93e-40 1.11e-75 0.00e+00 5.68e-63 1.00e+00 1.00e+00 1.00e-01
Optimal solution found
0.741599 seconds (1.53 M allocations: 118.708 MiB, 13.85% gc time)
iter time(s) μ P-obj D-obj gap P-error p-error d-error α_p α_d beta
Primal objective:119.9999999999999999999999999999999999999917627060555516439977164928930200937897
Dual objective:120.000000000000000000000000000000000000005990759232326077092569823350530853966
Duality gap:5.928355490322680456188887690629483406774762040292316850924073292798466630295741e-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.014130234s) **
** Transforming the problem and the solution ** (0.0031921980000000003s)
** Projection the solution into the affine space **
Reducing the system from 26 columns to 26 columns
Constructing the linear system... (0.099346616s)
Computing an approximate solution in the extension field... (0.004079824s)
Preprocessing to get an integer system... (0.000673588s)
Finding the pivots of A using RREF mod p... (0.001339153 0.001833129 s)
Solving the system of size 38 x 40 using the pseudoinverse... 0.002547738s
** Finished projection into affine space (0.112082842s) **
** Checking feasibility **
The slacks are satisfied (checked or ensured by solving the system)
Checking sdp constraints
done (0.0665268)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())
else
objective = Objective(0, Dict(), Dict())
push!(constraints, Constraint(obj-1, Dict(k => [1;;] for k=0:2d), Dict()))
endwhere 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.15831434739030571916973180897335312215428964315433748803634556059281377006215)