Error analysis

Math quiz

How much is

1 + 2 = \mathord{?}

1 + 2 = 3

print(1 + 2)  # 3

Math quiz

How much is

1.0 + 2.0 = \mathord{?}

1.0 + 2.0 = 3.0

print(1.0 + 2.0)  # 3.0

Math quiz

How much is

0.1 + 0.2 = \mathord{?}

0.1 + 0.2 = 0.3

print(0.1 + 0.2)  # 0.30000000000000004

Number representation

How do we represent numbers in a computer?

Computers only store bits, usually in groups of 8, 16, 32, or 64.

It is natural to represent integers in binary.

Loading diagram...

\underbrace{0} _ {\min} \le \underbrace{101101_2 = 23_{10}} _ {\text{value}} \le \underbrace{ 11\underset{27}{\cdots}11_2 = 2147483647_{10}} _ {\max}

Fixed-point

What about fractions?

A natural way is to use fixed-point numbers.

Loading diagram...

\begin{gather*} \underbrace{0} _ {\min} \le \underbrace{111001.101_2 = 23.625_{10}} _ {\text{value}} \le \\ \le \underbrace{ 11\underset{11}{\cdots}11.11\underset{12}{\cdots}11_2 = 32767.99998474121 _{10}} _ {\max} \end{gather*}

Pros-Cons of fixed-point numbers

Pros
- Stable error
- Simple operations
Cons
- Extremely Limited range
- Low precision

Floating-point

The IEEE 754 standard defines floating-point numbers.

Loading diagram...

\begin{gather*} \underbrace{0} _ {\min} \le \underbrace{1.101_2 \times 2^{2} = 6.5_{10}} _ {\text{value}} \le \\ \le \underbrace{ 1.11\underset{19}{\cdots}11_2 \times 2^{127} = (1 - 2^{-24}) _ {10} \times 2^{127}} _ {\max} \end{gather*}

Problems with floating-point

The number of bits is still finite.

Range: Limited exponent
Error: Rounding errors
Unit roundoff: $\epsilon$ is the smallest positive number such that $1 + \epsilon > 1$

Numbers with infinite decimals

Irrational numbers are not representable in a finite number of bits.

We always have to round them (approximation) or use symbolic computation.

Some fractions can create issues (think of $1/3$ ). They are called non-dyadic numbers.

Computation over floating point numbers

The standard dictates that the result of an operation is computed exactly and then rounded.

\begin{align*} a + b & \implies \textbf{fl}(a + b) = (a + b)(1 + \delta) \\ a - b & \implies \textbf{fl}(a - b) = (a - b)(1 + \delta)\\ a \times b & \implies \textbf{fl}(a \times b) = (a \times b)(1 + \delta) \\ a \div b & \implies \textbf{fl}(a \div b) = (a \div b)(1 + \delta) \end{align*}

with $|\delta| \le \epsilon$ .

Overflows and underflows are also possible, but they are detectable.

Errors in computation

Let the exponent and mantissa be 4 and 7 (+1) bits.

\underbrace{1.1111100_2 \times 2^{7}} _ {252_{10}} + \underbrace{1.1_2} _ {1.5_{10}} =\ ?

\begin{array}{crll} & 111111000 & + \\ & 11 & = \\ \hline & 111111011 & \\ \end{array}

The result would need 8 bits (+1) to be represented.

1.\underbrace{1111101} _ {\text{Stored}}1_2 \times 2^{7} = 253.0 \ne 252 + 1.5 = 253.5

Errors propagate

Errors propagate through operations.

\textbf{fl}(a + b + c) = ((a + b)(1 + \delta_1) + c)(1 + \delta_2)

The exact value depends on the order of operations, making floating point arithmetic non-associative.

Error analysis

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

Forward error analysis VS Backward error analysis.

Errors in inner products

The inner product of two column vectors is

\hat{s} = \textbf{fl}(x^Ty) = \textbf{fl}\left(\sum_{i=1}^{n} x_i y_i\right)

\hat{s} = x_1y_1(1 + \delta_{i\times})\prod^n_{i = 2}(i + \delta_i) + \sum_{i=2}^{n} x_i y_i(1 + \delta_{i\times})\prod^n_{j = i}(1 + \delta_i)

Typical notation

\begin{align*} &\prod^n_{i = 1}(1 + \delta_i)^{p_i} = 1 + \theta_n \\ &|\theta_n| \le \frac{n\epsilon}{1 - n\epsilon} \eqqcolon \gamma_n \end{align*}

Error bound

\textbf{fl}(x^Ty) \le x^Ty(1 + \gamma_n)

Using the backward error analysis,

\textbf{fl}(x^Ty) = (x + \Delta x)^Ty = x^T(y + \Delta y)

with

|\Delta x| \le \gamma_n|x|, \quad |\Delta y| \le \gamma_n|y|

LU decomposition

Compute the $LU$ factorisation of

A = \begin{bmatrix} \epsilon & -1 \\ 1 & 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ l_{21} & 1 \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} \\ 0 & u_{22} \end{bmatrix}

u_{11} = \epsilon, \quad u_{12} = -1, \quad l_{21} = \epsilon^{-1}, \quad u_{22} = 1 + \epsilon^{-1} \approx\epsilon^{-1}

\hat{L}\hat{U} = \begin{bmatrix} \epsilon & -1 \\ 1 & 0 \end{bmatrix}

Notable results

Triangular systems: $(T + \Delta T)\hat{x} = b \quad |\Delta T| \le \gamma_n|T|$
LU decomposition: $(A + \Delta A) = \hat{L}\hat{U} \quad |\Delta A| \le \gamma_n|\hat{L}||\hat{U}|$
LU system: $(A + \Delta A)\hat{x} = b \quad |\Delta A| \le \gamma_{3n}|\hat{L}||\hat{U}|$

Simplex algorithm

The simplex is the most widely used method for solving LP problems.

\begin{align*} \text{maximize} \quad & c^Tx \\ \text{subject to} \quad & Ax = b \\ & x \ge 0 \end{align*}

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

We define a basis $B$ of $m$ columns of $A$ . It must be non-singular.

Revised simplex

The revised simplex is the basis for the more common implementation.

Find an initial feasible solution (auxiliary LP).
Primal solution $x_B = B^{-1}b$ .
Dual solution $y_B = B^{-T}c_B$ .
Entering variable $e = \arg\max_{r_e > 0}(r_e \coloneqq c_e - A_e^Ty_B)$ .
Coefficient $d = B^{-1}A_e$ .
Leaving variable $l = \arg\min_{d_i > 0}(x_{B_i}/d_i)$ .
Update $B = \begin{bmatrix} B_1 & \cdots & B_{l-1} & B_{l+1} & \cdots A_e \end{bmatrix}$ .

Exit conditions

The steps of the simplex loop until

No initial feasible solution $\implies$ infeasible problem.
No entering variable $\implies$ optimal solution.
No leaving variable $\implies$ unbounded problem.

Bartels–Golub decomposition

Three linear systems involving $B$ each iteration. Why not use $LU$ decomposition?

Keep $L$ fixed and only update $U$ .

H^{(i)} = \begin{bmatrix} U^{(i-1)}_1 & \cdots & U^{(i-1)}_{l_{i-1} - 1} & U^{(i-1)}_{l_{i-1} + 1} & \cdots & U^{(i-1)}_m & L^{-1}A_e \end{bmatrix}

Bartels–Golub decomposition

We need to re-triangularise $H^{(i)}$ .

\begin{align*} B^{(i)} &= LH^{(i)} \\ &= LC^{(1)}C^{(2)}\cdots C^{(i)}\underbrace{C^{(i)-1}C^{(i-1)-1}\cdots C^{(1)-1}}_{\text{Gaussian elimination}}H^{(i)} = \\ &= L\underbrace{C^{(1)}C^{(2)}\cdots C^{(i)}}_{G^{(i)}}U^{(i)} = \\ &= LG^{(i)}U^{(i)} \\ \end{align*}

Stability

It can be proven that the Bartels–Golub decomposition is backward stable.

As long as the iterations do not grow too much, the error is bounded.

At some point it becomes more convenient to recompute the decomposition.

Exact floating point simplex

$\exists \tau, \epsilon \in \mathbb{R}$ , $x, \tau > 0$ , such that

def simplex(A, b, c, B_idx, tolerance, epsilon):
    m, n = A.shape
    for i in range(n!/(n-m)!):
        B = A[:, B_idx]
        L, G, U = bartels_golub(B)
        x_B = solve(L, G, U, b)
        y_B = solve(L.T, G.T, U.T, c[B_idx])
        r = c - A.T @ y_B
        if r <= tolerance: return x_B # Optimal
        e = argmax(r)
        d = solve(L, G, U, A[:, e])
        if d <= tolerance: return x_B # Unbounded
        l = argmin(x_B / d, where=d > tolerance)
        B_idx[l] = e

behaves exactly.

Takeaways

Working with integers
Increasing the precision and tolerance
Using rational arithmetic
Using symbolic computation

Error analysis

Math quiz

Math quiz

Math quiz

Number representation

Fixed-point

Pros-Cons of fixed-point numbers

Floating-point

Problems with floating-point

Numbers with infinite decimals

Computation over floating point numbers

Errors in computation

Errors propagate

Error analysis

Errors in inner products

Typical notation

Error bound

LU decomposition

Notable results

Simplex algorithm

Revised simplex

Exit conditions

Bartels–Golub decomposition

Bartels–Golub decomposition

Stability

Exact floating point simplex

Takeaways

Resources