delpi

Verified Linear Programming through Tolerance-Aware Precision Boosting

E. Casablanca

M. Sidaway

P. Zuliani

S. Soudjani

VSTTE 2026

Bits of computer architecture

Computers use bits to store data, including numbers.

Floating point
$0.2, -42.42, 100.10$

Rationals
$\frac{1}{2}, -\frac{3}{4}, \frac{3}{3}$

We only have a finite number of bits.

$\frac{1}{10} + \frac{2}{10} = \frac{3}{10}$

$0.1 + 0.2 = 0.3$

$\frac{1}{3} + \frac{2}{3} = 1$

$0.3333\ldots + 0.6666\ldots = 1.0?$

Let me quickly refresh your memory about some fundamentals in computer architecture. All data in a computer is stored as bits, but how they are interpreted and, as a result, how the CPU operates on them, can vary. When it comes to non-integer numbers, the most common representation is floating-point, which is directly supported by the underlying hardware. Hence, it is very efficient to operate upon. This said, nothing prevents us from implementing more complex data types, such as rational numbers. In that case, we are also responsible for implementing the operations on them, such as addition, multiplication, etc. But why does the distinction matter? The reason is that we only have a finite number of bits, and therefore we cannot represent rational numbers exactly, let alone irrational ones. To illustrate this, just take a look at the following example, using base 10 for simplicity. We may be lucky and get the same result with both representations, but, in general, we will never have the infinite bits needed to represent all values exactly, and therefore we have to settle for an approximation. -->

Motivation

Standard linear program (LP)

\min \{ c^\top x \mid Ax=b,\ x \ge 0 \}

with $A \in \mathbb{Q}^{m \times n}$ , $b \in \mathbb{Q}^m$ , $c \in \mathbb{Q}^n$ .

Linear Programming is a core component of:
- Optimization
- SMT Solving
- Program Analysis
Numerical errors undermine soundness.
- Commercial floating-point LP solvers can produce incorrect results.

Solutions

Rational arithmetic
- Easy to implement, exactness guarantee.
- Slow and memory-intensive.

Iterative refinement (SoPlex $^{1}$ $^{1}$ )
- Generally faster.
- Can struggle with numerical issues.

Precision boosting (QSopt_ex $^{2}$ $^{2}$ )
- Can solve difficult LPs.
- Limited theoretical justification.

[1] A. Gleixner, D. Steffy, K. Wolter (2016), Iterative Refinement for Linear Programming

[2] D. L. Applegate, W. Cook, S. Dash, D. G. Espinoza (2009), QSopt_ex Rational LP Solver

Goal

Loading diagram...

IEEE 754 Standard

Modern computers implement the IEEE 754 standard $^{3}$ for floating-point arithmetic.

\text{Floating-point number} = \pm m \cdot \beta ^{e-t}

where $\beta$ is the base, $t$ is the precision, $m$ is the integer significand and $e$ is the exponent.

Loading diagram...

[3] IEEE (2019) IEEE Standard for Floating-Point Arithmetic

Operations with floating-point numbers

The IEEE 754 standard provides the following guarantees for floating-point operations:

\textbf{fl}_{\epsilon}\left({x \odot y}\right) = (x \odot y)(1 + \mathbb{\varepsilon})^{\pm 1} \quad \text{with } |\mathbb{\varepsilon}| \le \epsilon .

where $\epsilon$ is the unit roundoff ( $2^{-53} \approx 1.1 \times 10^{-16}$ for double precision) and $\textbf{fl}_{\epsilon}\left({x \odot y}\right)$ is the floating-point computation of $x \odot y$ .

Any operation can introduce an error.

Error analysis

Error analysis is the study of how floating-point computations affect our results.

Forward error analysis VS Backward error analysis $^{4}$ .

[4] N. J. Higham (2002), Accuracy and Stability of Numerical Algorithms

Ensuring backward stability in the simplex

Given an initial sequence of indices $\mathcal{B}_1, \mathcal{B}_2, \dots, \mathcal{B}_m$ corresponding to the columns of $A$ , at each iteration:

Construct the basis: $B = \begin{bmatrix} A_{:,\mathcal{B}_1} & A_{:,\mathcal{B}_2} & \dots A_{:,\mathcal{B}_m} \end{bmatrix}$
Compute primal vector: $x_\mathcal{B} = B^{-1}b$
Compute dual vector: $y = B^{-T}c_\mathcal{B}$
Evaluate reduced costs: $r = c - A^\top y$
Determine entering variable.
Determine leaving variable: $d = B^{-1}A_{:,e}$ .
Pivot and update basis.

Bartels-Golub update

After the initial $LU$ decomposition of $B$ , the factor $U$ is updated $^{5}$ :

U^{(i)}

H^{(i+1)}

U^{(i+1)}

[5] R. H. Bartels (1971), A Stabilization of the Simplex Method

Error analysis

In the backward-error analysis formulation,

(B^{(i)} + \mathcal{E}^{(i)})x = b, \qquad x, b \in \mathbb{Q}^m,

we find that

\|\mathcal{E}^{(i)}\|_1 = O(\epsilon).

As long as the number of iterations is limited, the backward error is reasonably bounded.

Variable precision floating-point simplex

Loading diagram...

Characteristic constants

For a given input $A$ , $b$ , and $c$ , we can define the following constants:

$r^-_{\max}$ is the maximum negative across all reduced costs;
$d^+_{\min}$ is the minimum positive element across all $d$ s;
$u^{1\sim2}_{\min}$ is the smallest difference between the two smallest distinct values of $u = x_{\mathcal{B}}/d$ across all bases.

Computing them directly would require exploring all the possible bases.

Convergence of precision boosting

If $(1)$ and $(2)$ hold, then the floating-point simplex makes the same decisions as the exact rational simplex.

(1) \quad \tau < -r^-_{\max} / 2 \qquad \tau < d^+_{\min} / 2 \qquad \tau < u^{1\sim2}_{\min}

(2) \quad \|\varepsilon r\|_1 \le \tau \qquad \|\varepsilon d\|_1 \le \tau \qquad \|\varepsilon u\|_1 \le \min(\tau, u^{1\sim2}_{\min} - \tau) / 2

Delta-optimality

Theorem

Given $\delta > 0$ , a $\delta$ -complete LP algorithm returns exactly one of the following:

$\delta$ -optimal: iff it produces a $\delta$ -optimal certificate $(x_U, y_L)$ such that $0 \le c^\top x_U - b^\top y_L \le \delta$ ;
infeasible: iff the LP is infeasible;
unbounded: iff the LP is unbounded.

If $\delta = 0$ , then the algorithm is exact.

Precision boosting simplex

Loading diagram...

Delpi

Delta-complete Exact Linear ProgrammIng framework: Delpi.
C++ library with python interface.
QSopt_ex as a basis.
Implements the $\delta$ -complete procedure.
Can interface with SoPlex.

Results

Sloane–Stufken benchmarks $^6$ .
Dense and feasible.
Numerically challenging.
- Incorrect results from FP SOTA solvers.
Performance profile plot.

[6] N. J. A. Sloane, J. Stufken (1996), A linear programming bound for orthogonal arrays with mixed levels

Future Work

Integration into SMT solvers (work in progress!)
Stronger practical heuristics for choosing tolerances
Extensions to MILP

References

Thank you!