Rounding numerical solutions to exact optimal solutions

In certain situations there are reasons to believe that there is a 'nice' optimal solution. That is, the kernel of every optimal solution can be defined over a number field of low algebraic degree and low bit size. In that case, we can use the rounding procedure from [1] to obtain such a solution from a precise enough solution returned by a semidefinite programming solver. This is implemented in the function exact_solution.

See also the following examples:

ClusteredLowRankSolver.exact_solution — Function
exact_solution(problem::Problem, primalsol::PrimalSolution, dualsol::DualSolution; transformed=false, FF = QQ, g=1, settings::RoundingSettings=RoundingSettings(), monomial_bases=nothing)

Compute and return an exact solution to the problem, given a primal solution, dual solution and a field FF with approximate generator g. Return (success, exactdualsol) if transformed=false, and (success, pd_transformed_exactsolution, transformations) if transformed=true.

A short introduction to the rounding procedure

The rounding procedure consists of three steps. First a nice basis for the kernel of the provided solution is found; this defines the optimal face. In the second step, this basis is used as part of a basis transformation, to obtain a smaller semidefinite program where the affine hull of the constraints has a one-to-one correspondence with the original optimal face. The provided solution is also transformed, to obtain a strictly feasible solution of the new semidefinite program. In the last step, we take an approximation of the numerical, transformed solution over the rationals (or a provided number field), and find an exact solution close to it. See [1] for a more detailed description.

Settings for the rounding procedure

For ease of use, all settings which can be used to tweak the rounding procedure are collected in the RoundingSettings structure.

ClusteredLowRankSolver.RoundingSettings — Type
RoundingSettings(settings...)

Settings for the rounding procedure:

  1. Finding the kernel
    • kernel_errbound: (default: 1e-10) the allowed error for the kernel vectors. That is, the maximum entry of Xv is in absolute value at most this
    • kernel_round_errbound: (default: 1e-15) the maximum allowed error when rounding the reduced row-echelon form entrywise
    • kernel_use_primal: (default: true) use the primal solution to find the kernel vectors (otherwise, use an SVD of the dual solution)
    • kernel_lll: (default: false) if true, use the LLL algorithm to find the nullspace of the kernel. Otherwise, use the reduced row-echelon form to find kernel vectors.
    • kernel_bits: (default: 1000) the maximum number of bits to be used in the LLL algorithm (for finding relations or finding the entries of the RREF)
  2. Reducing the kernel vectors
    • reduce_kernelvectors: (default: true) apply the reduction step or not
    • reduce_kernelvectors_cutoff: (default: 400) do reduction on the full matrix if its size is at most this cutoff. Otherwise do it on a submatrix
    • reduce_kernelvectors_stepsize: (default: 200) the number of extra columns to take into account in each iteration of the reduction step
  3. Transforming the problem and the solution
    • unimodular_transform: (default: true) use a unimodular transform obtained in the reduction step
    • normalize_transformation: (default: true) multiply by a diagonal matrix to get an integer transformation for the problem (for problems over QQ)
  4. Finding an approximate solution in the field
    • regularization: (default: 1e-20) use this regularization for solving the extended system
    • approximation_decimals: (default: 40) Approximate the numerical solution using this many digits, entrywise
  5. Rounding the solution to the affine space of constraints
    • redundancyfactor: (default: 10) take at least this times the number of constraints columns as potential pivots
    • pseudo: (default: true) use the psuedo inverse for rounding (this may give solutions with larger bitsize than needed)
    • pseudo_columnfactor: (default: 1.05) For a system of r rows, use r * pseudo_columnfactor number of columns for the pseudo inverse
    • extracolumns_linindep: (default: false) if true, take the extra columns linearly independent of each other (otherwise, random columns)

Finding the right settings

Depending on the problem, the default parameters will suffice. In the following cases small changes are needed:

  1. The kernel vectors are not found correctly, or have significantly higher maximum number after reduction than the maximum numerator and denominator.
    • Try to solve to a smaller duality gap, and/or decrease the setting kernel_round_errbound.
    • Try the settings kernel_use_primal=false.
  2. The final solution does not satisfy the constraints. (Or not enough pivots were found)
    • increase redundancyfactor, or set it to -1 to take all variables into account.
    • This might also be an indication that the kernel vectors are not found correctly (see item 1).
  3. The final solution is not positive semidefinite.
    • Increase the setting approximation_decimals. Make sure that the provided solution has at least that many digits correct (the duality gap should be lower than 10^(-approximation_digits))
    • Increase the setting pseudo_columnfactor. The higher this setting, the closer the exact solution will be to the provided solution. However, this also increases the bit size of the exact solution.
    • In some cases this is also an indication that the kernel vectors are not found correctly (see item 1), especially when the reported negative eigenvalues are close to zero.

Using coefficient matching

Although the semidefinite program used in the solver is defined using sampling, it is possible to use a semidefinite program defined using coefficient matching in a monomial basis for the rounding procedure. This is not yet fully automated; the user needs to provide the monomial basis for each polynomial constraint to the rounding procedure in order to use this. Using coefficient matching instead of sampling in the rounding heuristic generally decreases the size of the exact solutions. See below for an example using the slightly modified code for the Delsarte LP bound in the example for rounding. Here we have one univariate constraint with polynomials up to degree $2d$.

d = 3
problem, primalsol, dualsol = delsarte(8, 1//2, d; duality_gap_threshold=1e-30)
R, (x,) = polynomial_ring(QQ, 1)
mon_basis = [x^k for k=0:2d]
success, exactdualsol = exact_solution(problem, primalsol, dualsol, monbases = [mon_basis])

Finding the appropriate number field for the rounding procedure

It is in general not directly clear over which number field the optimal face can be defined. The function find_field can help to find such a field. See Section 2.5 of [1] for an explanation of the procedure.

ClusteredLowRankSolver.find_field — Function
find_field(primalsol, dualsol, max_degree=4; valbound=1e-15, errbound=1e-15, bits=100, max_coeff=1000)

Heuristically find a field over which the kernel can probably be defined.

Only consider values at least valbound in absolute value. Find minimal polynomials such that the chosen entries are approximately generators with an error bound of errbound. Use bits number of bits and reject minimal polynomials with a maximum coefficient of more than max_coeff.

Getting an exact feasible solution close to the optimum

In general, the kernel of the optimal solutions can only be defined using a number field of high algebraic degree, or with large bit size. In this case it is unlikely that the exact kernel can be found, so that the rounding procedure will fail. However, it is possible to get an exact feasible solution with objective value close to the optimum using the rounding procedure. To do this, solve the semidefinite program as a feasibility problem (no objective), with the extra constraint that the (original) objective should be a rational number close to the optimum value. Then the solver will return a strictly feasible solution if it exists (that is, with empty kernel). The function exact_solution will now essentially skip the steps of finding the kernel vectors and transforming the problem, since there are no kernel vectors. Therefore, it will find an exact feasible solution close to the provided solution such that the objective of the original problem is equal to the rational number in the extra constraint. See Rounding the Delsarte LP bound for a worked example.