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

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$
Infinite digits: Irrational numbers or non-dyadic fractions

Below the unit roundoff

#include <stdio.h>

int main() {
    float meters = 0;
    int iterations = 100000000;
    for (int i = 0; i < iterations; i++) {
        meters += 0.01;
    }
    printf("Expected: %f km\n", 0.01 * iterations / 1000 );
    printf("Got: %f km \n", meters / 1000);
}

Output:

// Expected: 1000.000000 km
// Got: 262.144012 km

The gap is too large: $[262144.0], [262144.03125], [262144.0625]$

The importance of partial pivoting

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}

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} + \underbrace{1.1_2} _ {1.5} =\ ?

\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 in inner products

The inner product of two column vectors is

x, y \in \mathbb{R}^n \newline \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^{\times}_{i})\prod^n_{i = 2}(i + \delta^{+}_i) + \sum_{i=2}^{n} x_i y_i(1 + \delta^{\times}_{i})\prod^n_{j = i}(1 + \delta^{+}_i)

Error bound

\hat{s} = \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|

Notable results

Triangular systems: $(T + \Delta T)\hat{x} = b \quad |\Delta T| \le \gamma_n|T|$
LU decomposition (PP): $(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 among the most widely used methods 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}$

Bartels–Golub decomposition

Three linear systems involving $B$ each iteration.

$x_B = B^{-1}b$ .
$y_B = B^{-T}c_B$ .
$d = B^{-1}A_e$ .

Use $LU$ decomposition. Keep $L$ fixed and only update $U$ .

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

Bartels–Golub decomposition

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

\begin{align*} B^{(i)} &= LG^{(i - 1)}H^{(i)} \\ &= LC^{(1)}C^{(2)}\cdots C^{(i)}\underbrace{\Gamma^{(i)}_{m-1}{\Pi}^{(i)}_{m-1}\dots {\Gamma}^{(i)}_{l_i}{\Pi}^{(i)}_{l_i}}_{C^{(i)-1}}H^{(i)} = \\ &= L\underbrace{C^{(1)}C^{(2)}\cdots C^{(i)}}_{G^{(i)}}U^{(i)} = \\ &= LG^{(i)}U^{(i)} \\ \end{align*}

with $i > 0, \quad G^{(0)} \coloneqq I$ .

Error analysis

Error analysis

Math quiz

Math quiz

Math quiz

Number representation

Fixed-point

Pros-Cons of fixed-point numbers

Floating-point

Problems with floating-point

Below the unit roundoff

The importance of partial pivoting

Computation over floating point numbers

Errors in computation

Errors propagate

Error analysis

Errors in inner products

Typical notation

Error bound

Notable results

Simplex algorithm

Revised simplex

Exit conditions

Bartels–Golub decomposition

Matrix shapes

Bartels–Golub decomposition

Stability

Exact floating point simplex

Exact rational simplex

Key steps

Advantages

Takeaways

Resources