# What approximation is guarantees when solving an LP with floating-point numbers?

Given a linear program

\begin{align} \text{maximize} \quad & c^{T}x \\ \text{s.t.} \quad & A x \leq b \end{align}

I can solve it exactly in polynomial time, using e.g. interior-point methods.

But in practice, most solvers perform all computations using floating-point numbers. For concreteness, suppose I use floating-point numbers with 32 bits. What guarantees can be formally proved on the returned solution?

For example, I would be happy to have a theorem such as "for every solution $$x$$ returned by the solver:

• The constraints are approximately satisfied, i.e., $$A x \leq b + 10^{-10}$$;
• The solution is approximately optimal, i.e., for every $$y$$ that satisfies the constraints, $$c x + 10^{-20} \geq c y$$; ".

Does there exist such a theorem (with possibly different constant values), or a similar one?

• I would be surprised if you could find such a guarantee, since it is known that "solutions" can be horribly wrong if a basis matrix has an ugly condition number.
– prubin
Commented Apr 10 at 18:18

A quick disclaimer regarding

I can solve it exactly in polynomial time using e.g. interior-point methods

• Interior-point algorithms do not solve a problem "exactly": they provide a solution up to a $$\epsilon > 0$$ accuracy. Note that this is an algorithmic property, i.e., it has nothing to do with finite-precision arithmetic. Obtaining an exact solution to the LP problem would require something like a crossover step.
• Interior-point methods provide an $$\epsilon$$-optimal solution in $$O(\log 1/\epsilon)$$, but they are not strongly polynomial-time for solving an LP to exact optimality (in addition to the limitation above). See, e.g., this paper and this paper.
• When computing in finite-precision arithmetic, as far as I know, only LPs with rational input data (i.e., all coefficients of $$A, b, c$$) can be solved exactly. Typically this is done using a combination of simplex, rational arithmetic, and possibly iterative refinement (see, e.g., this recent paper for a nice overview).

For simplicity, I will consider equality constraints $$Ax = b$$ in what follows. Similar statements apply to inequality constraints.

That being said, assuming you are using an existing LP solvers with finite-precision arithmetic (64-bit floating-point precision in 99.99% of cases nowadays). It's important to note that all solvers have numerical tolerances. Assuming no numerical issue, the returned solution is only "feasible" up to the solver's tolerances.

• Simplex solvers typically use an absolute tolerance $$\epsilon_{abs}$$ of about $$10^{-6}$$, i.e., $$x$$ is considered feasible if $$\|Ax - b\|_{\infty} \leq \epsilon_{abs}$$.
• Interior-point solvers typically use a relative tolerance $$\epsilon_{rel}$$ of about $$10^{-8}$$, i.e., for an equality constraint $$Ax = b$$, $$x$$ is considered feasible if $$\|Ax - b\|_{\infty} \leq \epsilon_{rel} \|b\|_{\infty}$$.

Note that the numerical tolerances I give here are not set in stone, they are the default values for solvers I know, e.g., CPLEX, Gurobi, Mosek, etc. Also note that tolerances are typically enforced on the scaled problem, which may yield larger violations on the original problem passed by the user.

As far as I know, no solver is guaranteed to always return a solution that satisfies user-prescribed tolerances, because of numerical issues. The solver may return a solution that satisfies the tolerances, it may return a reduced-accuracy solution, or no solution at all. Furthermore, as mentioned by prubin in his comment, small residuals do not imply much regarding 1) distance to feasibility and 2) distance to optimality.

Namely, $$\|Ax - b\|$$ being small does not imply that the distance between $$x$$ and the feasible set is small (in absolute terms), nor does it imply that $$x$$ is close to optimality. For instance, consider the problem $$\min_{x} \{ x \ | \ \rho x = \rho \}$$, for $$\rho > 0$$. There is only one feasible (and therefore optimal) solution: $$x = 1$$. However, if $$\rho = 10^{-20}$$, then $$x = 10^6$$ has residual roughly $$10^{-14}$$, which may be considered "small" in absolute standards. However, $$10^{6}$$ is obviously far from optimum. This small example illustrates prubin's comment that the condition number of the constraint matrix $$A$$ can have a dramatic influence on solution quality, even in the presence of "small" residuals.

• Is there any relation between the absolute/relative tolerance of the solver, and the length of floating-point numbers that I use? E.g. if I use 128-bit instead of 64-bit, will the tolerance be smaller? Commented Apr 11 at 18:19
• It will depend on what tolerances you set in the solver, assuming said solver supports higher-precision arithmetic. My rule of thumb is to use the square root of the numerical precision of the arithmetic. That's roughly 1e-8 for 64-bit, and 1e-16 for 128-bit. Those are the default tolerances in Tulip: ds4dm.github.io/Tulip.jl/stable/reference/options/#Tolerances. Disclaimer: I wrote that solver and it hasn't been updated in a while. Commented Apr 11 at 23:30