Reliability of ILP approach to number-theoretic optimization

Question

This question is inspired by the recent answer, where @RobPratt proposed to use integer linear programming (ILP) for solving a number-theoretic optimization problem. I will consider a very similar problem described by sequence OEIS A061057, which for a given integer $n$ asks to find a split of $n!$ into two co-factors with the smallest difference, i.e. to find $$\min_{d\,\mid\, n!}\quad \left|\frac{n!}{d} - d\right|.$$ Equivalently, it can be posed as finding the largest divisor $d\mid n!$ such that $d\leq \sqrt{n!}$.

---

This problem has a simple ILP formulation based on the prime factorization $$n! = p_1^{m_1} p_2^{m_2} \cdots p_k^{m_k}$$ with integer variables $x_1,\dots,x_k$, representing the exponents in the prime factorization of $d=p_1^{x_1}p_2^{x_2}\cdots p_k^{x_k}$, as follows: $$\begin{cases} \sum_{i=1}^k x_i\log p_i\ \longrightarrow\ \max,\\ \sum_{i=1}^k x_i\log p_i\leq \log \sqrt{n!},\\ 0\leq x_i\leq m_i,\quad i\in\{1,2,\dots,k\}. \end{cases}$$

---

Existing ILP solvers are pretty happy with this problem formulation and solve it for $n\leq 40$ in a matter of seconds. The issue, however, is that their solutions quickly become incorrect.

I've tried to solve this problem with 3 solvers: Gurobi, CPLEX, and GLPK, with default parameters from Sage. The code using GLPK can be even run online at SageMathCell. As an example, it computes the difference (= 928) for $n=15$.

I've compared their results to the known values from A061057. To my surprise, the celebrated commercial solvers Gurobi and CPLEX showed worst and quite similar performance by producing incorrect results for $n\geq 18$ and $n\geq 19$ respectively. The freeware GLPK did much better job by failing only for $n=27$ and $n\geq 31$.

Gurobi has an extensive discussion on numerical issues, which does not seem to be much relevant here, especially given that GLPK's performance was significantly better.

---

So, I have a few related questions:

Q1. What is the cause of failure for ILP solvers in this problem?

Q2. Is there a way to alleviate the issue and extend the range of correctly soluble $n$?

Q3. Assuming that the correct answer is unknown, how must trust we can put in the result produced by an ILP solver and/or how we can verify it?

Answer 1

Well, the question is broad, but let us address at least some of it.

Q1. What is the cause of failure for ILP solvers in this problem?

Let's concentrate in the GLPK failure with $n=27$. Here is a SageMath cell for demonstrating what happens, based on your original code with some extra outputs. And here is the output.

we were using absolute tolerance 1e-05
we were using relative tolerance 1e-07
objective should <=  32.278769313503155
objective is         32.27878369954806
absolute exceedance  1.4386044902892081e-05
relative exceedance  4.4568133199780874e-07
(-3002360256, {2: 6, 3: 6, 5: 6, 7: 1, 11: 2, 13: 2}, {2: 17, 3: 7, 7: 2, 17: 1, 19: 1, 23: 1})

We note that the offered solution is invalid; the first factor is $2^6 3^6 5^6 7^1 11^2 13^2 = 104351247000000$, which is more than $\sqrt{27!}$.

Numerically, we were trying to maximize the objective value $\sum_i x_i \log p_i$, while simultaneously constraining it to be at most $\log \sqrt{27!}$, represented as a float. Observe that the resulting objective exceeds the constraint. The solver seems to be happy because the constraint is almost satisfied. (It seems like it is slightly outside the tolerance, but it is pretty near, and I did not read enough of the documentation to be 100% sure that these are the only tolerance parameters affecting the behaviour.)

I think this is an instructive example: The failure is not caused by rounding errors when converting the original inputs into floats (although they could be a problem too). It is caused by the relatively loose-looking tolerance in the solver. Such tolerance values are in fact quite common in numerical solvers. And they might come as a big surprise if you are expecting exact results; especially since the user did not ask for such tolerance, it was just the solver's default.

You might hope that you could simply reduce the tolerance. But the GLPK manual warns: "(Do not change this parameter without detailed understanding its purpose.)"

---

Q2. Is there a way to alleviate the issue and extend the range of correctly soluble n?

I would try either of the following options:

1. Use stricter tolerance and higher precision floats, if the solver allows this. Then read the solver's manual very carefully, to find what (if anything) it actually guarantees for the solutions. Then double-check the solutions (e.g. are they even feasible, in the exact mathematical formulation of the problem). But even if we rigorously validate that the solver's numerical solution is feasible, do we really know it was optimal? Well, we can hope...

2. Use an exact solver (like PPL in my comments). Since the problem involves logarithms of rational numbers, we must approximate them with rationals, within some precision. How much precision to needed, is patently problem-specific, but at least in this problem it should be doable. Fix some precision and work out the worst possible error that the rounding could cause, and if this is small enough not to change the integer results, we are happy. At least if we have confidence on the solver...

Detailed solution for this problem

In this particular problem (factorial-splitting), there are only two places where we must approximate reals with rationals: the upper bound $\log\sqrt{n!}$ and the logarithms of the prime factors. We can simply make sure that the former is approximated upwards and the latter are approximated downwards. Then we will find a rational solution with objective value at least as big as the true optimal value. We might find something where the objective value in fact exceeds the true upper bound, but we can check that afterwards. If that happens, increase the precision until we find a feasible solution. Then we know we got the true optimum (assuming, of course, that we trust the computer and the ILP solver).

Here is a SageMath cell demonstrating the exact rational solution with $20!$, and what happens if the precision is not enough: we get an error message instead of an incorrect solution. So if a solution is found, I would be fairly confident it is correct. Here are the outputs with 10, 20 and 30 bits precision:

using 10-bit approximation
done in 1.243 s
Error: First factor too big
None

using 20-bit approximation
done in 1.042 s
(20, 800640, {2: 7, 3: 2, 5: 3, 7: 2, 13: 1, 17: 1}, {2: 11, 3: 6, 5: 1, 11: 1, 19: 1})

using 30-bit approximation
done in 1.061 s
(20, 800640, {2: 7, 3: 2, 5: 3, 7: 2, 13: 1, 17: 1}, {2: 11, 3: 6, 5: 1, 11: 1, 19: 1})

And here are the results for $n=40,\ldots,45$. OEIS has them up to $n=41$ and they match.

(40, 470500040794291200, {2: 10, 3: 15, 5: 2, 7: 1, 11: 1, 13: 3, 19: 1, 23: 1, 29: 1, 31: 1, 37: 1}, {2: 28, 3: 3, 5: 7, 7: 4, 11: 2, 17: 2, 19: 1})
(41, 2323929740464193400, {2: 35, 3: 3, 5: 2, 7: 4, 11: 1, 17: 1, 19: 1, 23: 1, 31: 1, 41: 1}, {2: 3, 3: 15, 5: 7, 7: 1, 11: 2, 13: 3, 17: 1, 19: 1, 29: 1, 37: 1})
(42, 20720967220237197312, {2: 26, 3: 13, 7: 3, 13: 3, 17: 1, 23: 1, 29: 1, 41: 1}, {2: 13, 3: 6, 5: 9, 7: 3, 11: 3, 17: 1, 19: 2, 31: 1, 37: 1})
(43, 69638496398882611200, {2: 13, 3: 10, 5: 7, 7: 3, 11: 1, 13: 1, 17: 2, 19: 2, 31: 1, 41: 1}, {2: 26, 3: 9, 5: 2, 7: 3, 11: 2, 13: 2, 23: 1, 29: 1, 37: 1, 43: 1})
(44, 61690805562507264000, {2: 20, 3: 5, 5: 6, 7: 6, 11: 2, 13: 2, 17: 2, 19: 1, 31: 1}, {2: 21, 3: 14, 5: 3, 11: 2, 13: 1, 19: 1, 23: 1, 29: 1, 37: 1, 41: 1, 43: 1})
(45, 416497216789463040000, {2: 24, 3: 2, 5: 6, 7: 3, 11: 1, 13: 3, 17: 1, 23: 1, 29: 1, 31: 1, 37: 1, 43: 1}, {2: 17, 3: 19, 5: 4, 7: 3, 11: 3, 17: 1, 19: 2, 41: 1})

---

Q3. Assuming that the correct answer is unknown, how must trust we can put in the result produced by an ILP solver and/or how we can verify it?

If the solver is exact, and if you have verified that your input encoding does not cause spurious results (see above), the remaining suspects are the solver itself (both algorithm and implementation), and hardware (random bit flips in the memory, design errors like the infamous Pentium FDIV bug). Some of these suspicions could be alleviated by independent runs (different hardware, different solver, different parameters).

Another prospect is certificates of optimality. You could hope that the ILP solver could output such a certificate and it could be independently verified. I'm not aware of any common ILP solver that does this, but there is some research towards this direction. Here are just a few hits from googling:

Cheung, Kevin K. H.; Moazzez, Babak, Certificates of optimality for mixed integer linear programming using generalized subadditive generator functions, Adv. Oper. Res. 2016, Article ID 5017369, 11 p. (2016). ZBL1387.90153.

Cheung, Kevin K. H.; Gleixner, Ambros; Steffy, Daniel E., Verifying integer programming results, Eisenbrand, Friedrich (ed.) et al., Integer programming and combinatorial optimization. 19th international conference, IPCO 2017, Waterloo, ON, Canada, June 26–28, 2017. Proceedings. Cham: Springer. Lect. Notes Comput. Sci. 10328, 148-160 (2017). ZBL1418.90176.

Answer 2

I would suggest using exact ILP solvers to rule out errors that are rooted in limited precision.

QSopt_ex for example works with exact rational numbers.

Another more recent exact solver is SCIP of the Zuse Institute in Berlin.