delpi

Floating point
$0.2, -42.42, 100.10$

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

$\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?$

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}$)
  • Generally faster.
  • Can struggle with numerical issues.

• Precision boosting (QSopt_ex$^{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:

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

Bartels-Golub update

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

[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

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

Precision boosting simplex

Loading diagram...

Delpi

Results

[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!