Modern error analysis in finite-dimensional settings traces its origins to the seminal work of James Hardy Wilkinson. It distinguishes between forward error, which is the error in the result of an algorithm, and backward error, which measures the perturbation in the input data.
In the following, we denote by
:label: ex-error-analysis-taylor
Let $f : \mathbb{R}^n \to \mathbb{R}$ be a differentiable function and
$$
\tilde{f}(x) = f(\bar{x}) + \nabla f(\bar{x}) \cdot (x - \bar{x})
$$
its first-order Taylor expansion about $\bar{x}$. We use the convention here that $\nabla f$ is a row and $x - \bar{x}$ a column vector.
:label: ex-error-analysis-difference-quotient
Approximate the derivative $f := g'$ by a difference quotient of $g$:
$$
\tilde{f}(x) = \frac{g(x + h) - g(x)}{h}.
$$
for some $h > 0$.
:label: ex-error-analysis-linear-system
Consider the solution of the non-singular linear system of equations
$$
Ay = b
$$
with $A \in \mathbb{R}^{n \times n}$ and $b \in \mathbb{R}^n$. Then, the input space is the product space $D = X = \mathbb{R}^{n \times n} \times \mathbb{R}^n$ and our input data is the tuple $(A, b) \in X$. The output space is $Y = \mathbb{R}^n$. The function $f$ is given as
$$
f(A, b) = A^{-1}b
$$
and the approximate function $\tilde{f}$ is, for example, the Gaussian elimination algorithm with rounding errors to compute the solution $y$.
Let
The relative forward error
whenever
The backward error addresses the following question:
Given an input
$x$ and the corresponding output$\tilde{y} = \tilde{f}(x)$ from the perturbed function$\tilde{f}$ , by how much would one need to perturb$x$ to an input$\tilde{x} = x + \Delta x$ so that the original function returns$\tilde{y}$ ? Thus we look for
$$ f(\tilde{x}) = \tilde{y} = \tilde{f}(x). $$
In other words, we ask if the exact function
:label: ex-backward-error-linearisation
Let $f(x) = x^2 + x$ and $\tilde{f} = x$, its linearisation at $0$. Let $x = 1$ so that $\tilde{y} = 1$. The function $f$ attains $1$ at $\tilde{x}_{1,2} := - \frac{1}{2} \pm \frac{\sqrt{5}}{2}$.
```{image} ./figures/backward_error.png
:width: 400px
:align: center
```
Therefore, $\Delta x_1 := \tilde{x}_1 - x = (- \frac{1}{2} + \frac{\sqrt{5}}{2}) - 1$ and $\Delta x_2 := \tilde{x}_2 - x = (- \frac{1}{2} - \frac{\sqrt{5}}{2}) - 1$ are perturbations of $x$ to reproduce $\tilde{f}$'s output.
:label: ex-backward-error-reversed-roles
In the second example, we reverse the roles of $f$ and $\tilde{f}$. Let $f(x) = x$ and $\tilde{f} = x^2 + x$. Given $\tilde{y} = 1$, there are two possible inputs: $- \frac{1}{2} \pm \frac{\sqrt{5}}{2}$. Re-reading the question in the box above clarifies that we assume the choice of input is given: suppose it is the positive input $x = - \frac{1}{2} + \frac{\sqrt{5}}{2}$. Now, with $\tilde{x} = 1$, the perturbation is $\Delta x_1 := \tilde{x} - x = 1 - (- \frac{1}{2} + \frac{\sqrt{5}}{2})$.
The potential lack of bijectivity of
Therefore, the absolute backward error of
To express the above, we can also sometimes write more compactly
Similarly, the relative backward error of
Recall that
The definitions simplify if we assume that both
Analogously, the minimum of
Thus, for invertible
:label: ex-backward-error-linear-systems-noninvertible
Let $f(A, b) = A^{-1} b$ as introduced in this section. Then, given a vector $\tilde{y}$, $f^{-1}(\tilde{y})$ is the set of matrix-vector pairs $(A,b)$ such that $A \tilde{y} = b$. Clearly, there are many linear systems that have the solution $\tilde{y}$ and therefore this $f$ is not invertible.
The definition of the backward error may seem convoluted at first, but it is very useful when combined with the condition number.
We assume momentarily that
Equality (2) is possible by choosing
Inequality (5) demonstrates that the absolute forward error depends on a product of two quantities: the sensitivity of
The absolute condition number of
with
Similarly, the relative condition number is defined as
while the local relative condition number is, for
Hence,
We use here gradient and Jacobian as interchangeable concepts. Returning to equation
for
and, for
In summary, the relative forward error is bounded by the product of the condition number and the relative backward error.
What is a good condition number in practice? Ideally, we want the condition number to be as small as possible so that errors in the input data are not amplified in the output error. For each magnitude of the condition number, we lose one digit of accuracy in the output error. This can be seen as follows: Assume we work in double precision and that the algorithm has a backward error that is a small multiple of
If a condition number is small, we say that the problem is well-conditioned. If the condition number is large, we say that the problem is ill-conditioned. The precise meaning of these terms depends on the requirements of the application.
The condition number is a worst-case bound. In practice, the actual forward error may be better than predicted by the condition number.
In the following, we will generally not use the subscript
When working with numerical computations, derivatives are frequently required. The above formulas for the condition numbers are one example. There are several ways to approximate or compute them. We begin with a numerical approximation by finite differences, followed by alternative approaches using symbolic and automatic differentiation.
An efficient way to approximate the derivative of a function
Choosing a suitable step size
import numpy as np
def f(x):
return np.sin(x**2)
def finite_difference_derivative(f, x, h=1e-5):
return (f(x + h) - f(x - h)) / (2*h)
# Approximate f'(1)
approx_derivative = finite_difference_derivative(f, 1.0)
print("Finite difference approximation of f'(1):", approx_derivative)Python also has libraries, such as sympy, that can compute derivatives symbolically, similar to manual calculations or computer algebra systems such as Mathematica.
import sympy as sp
# Define a symbolic variable
x = sp.Symbol('x', real=True)
# Define the function f(x) = sin(x^2)
f = sp.sin(x**2)
# Compute the derivative f'(x)
f_prime = sp.diff(f, x)
# Print the result
print("f'(x) =", f_prime)
# Optionally, evaluate the derivative at a specific point, e.g., x = 1
print("f'(1) =", f_prime.subs(x, 1))Another way to compute derivatives in Python is through libraries that implement automatic differentiation at runtime, rather than relying on symbolic manipulation. One such tool is autograd, which can take ordinary Python functions using NumPy and return functions representing their derivatives. This approach is particularly convenient for more complex functions or situations where symbolic differentiation is cumbersome.
import autograd.numpy as np
from autograd import grad
# Define the function f(x) = sin(x²)
def f(x):
return np.sin(x**2)
# Use autograd's grad to create a function for f'(x)
f_prime = grad(f)
# Evaluate the derivative at a specific point, for example, x = 1.0
print("f'(1) =", f_prime(1.0))The example shows how quickly one can obtain numeric derivatives without explicitly working out the derivative formula. The grad function tracks operations performed on x and automatically applies the chain rule, enabling fast and accurate derivative calculations.
:class: tip
Let us consider the function $f : D \subset \mathbb{R} \to Y, x \mapsto \sqrt{x}$ with $Y = \mathbb{R}$. Let $\| \cdot \|_X = \| \cdot \|_Y = | \cdot |$, i.e. the norms are equal to the modulus.
1. Show that $K_{abs} = 1/2$ if $D = [1,2]$.
2. Show that $K_{abs} = \infty$ if $D = [0,1]$.
3. Show that $\kappa_{rel} = 1/2$ if $D = (0,\infty)$.
**Part 1.** The absolute condition number is given by:
$$
K_{abs} = \sup_{x, x+\Delta x \in D \atop 0 < \| \Delta x \|} \frac{|f(x+\Delta x) - f(x)|}{|\Delta x|}.
$$
Here $f(x) = \sqrt{x}$ and $D = [1, 2]$. Using the mean value theorem, there exists $c \in (x, x+\Delta x)$ such that:
$$
f(x+\Delta x) - f(x) = f'(c) \Delta x,
$$
where $f'(x) = \frac{1}{2\sqrt{x}}$. Substituting:
$$
\frac{|f(x+\Delta x) - f(x)|}{|\Delta x|} = |f'(c)| = \frac{1}{2\sqrt{c}}.
$$
To find the supremum of $\frac{1}{2\sqrt{c}}$ for $c \in [1, 2]$, observe that $\frac{1}{2\sqrt{c}}$ decreases as $c$ increases. The maximum occurs at $c = 1$:
$$
K_{abs} = \frac{1}{2\sqrt{1}} = \frac{1}{2}.
$$
**Part 2.** For $D = [0, 1]$, consider:
$$
\frac{|f(x+\Delta x) - f(x)|}{|\Delta x|} = \frac{1}{2\sqrt{c}},
$$
where $c \in (x, x+\Delta x)$ and $c \in [0, 1]$. As $c \to 0^+$, $\sqrt{c} \to 0$, and $\frac{1}{2\sqrt{c}} \to \infty$. Thus, $K_{abs} = \infty$.
**Part 3.** The relative local condition number is:
$$
\kappa_{rel}(x) = \lim_{\delta \to 0} \sup_{x, x+\Delta x \in D \atop 0 < |\Delta x| \leq \delta} \frac{|f(x+\Delta x) - f(x)| / |f(x)|}{|\Delta x| / |x|}.
$$
Simplify:
$$
\frac{|f(x+\Delta x) - f(x)| / |f(x)|}{|\Delta x| / |x|} = \frac{|f(x+\Delta x) - f(x)|}{|\Delta x|} \cdot \frac{|x|}{|f(x)|}.
$$
From earlier, $\frac{|f(x+\Delta x) - f(x)|}{|\Delta x|} = \frac{1}{2\sqrt{c}}$, and $f(x) = \sqrt{x}$, so $\frac{|x|}{|f(x)|} = \sqrt{x}$. Substituting:
$$
\frac{|f(x+\Delta x) - f(x)| / |f(x)|}{|\Delta x| / |x|} = \frac{1}{2\sqrt{c}} \cdot \sqrt{x}.
$$
As $\delta \to 0$, $c \to x$, so:
$$
\kappa_{rel}(x) = \frac{\sqrt{x}}{2\sqrt{x}} = \frac{1}{2}.
$$
:class: tip
Assuming single precision floating-point arithmetic, consider an algorithm with a backward error of $2 \cdot \epsilon_{mach}$ for a problem with a condition number $\kappa = 10^{3}$. How many correct digits do you expect in your solution?
For single precision, we have $\epsilon_{mach} \approx 6 \times 10^{-8}$. The condition number is $10^3$. Hence, the forward error is bounded by $1.2 \times 10^{-4}$, giving us nearly 4 digits of accuracy in the solution.
:class: tip
Let $f : D \subset X \to Y$ be differentiable, where $X$ and $Y$ are finite-dimensional real vector spaces. Show that
$$
\kappa_{abs}(x) = \| \nabla f(x) \|_{Y,X}, \qquad \kappa_{rel}(x) = \frac{\| x \|_X}{\| f(x) \|_Y}\| \nabla f(x) \|_{Y,X}.
$$
Hint: You may find it easier to first try the case $X = Y = \mathbb{R}$ with $\| \cdot \|_X = \| \cdot \|_Y = \| \cdot \|_{Y,X} = | \cdot |$.
Because $X$ and $Y$ are finite-dimensional real vector spaces, we may identify them as $X = \mathbb{R}^n$ and $Y = \mathbb{R}^m$. Noting that $f$ is differentiable at $x$, there is a function $h: \mathbb{R}^n \to \mathbb{R}^{m \times n}$ with $\lim_{\| \Delta x \| \to 0} h(x + \Delta x) = 0$ such that
$$
f(x + \Delta x) = f(x) + \nabla f(x) \cdot \Delta x + h(x + \Delta x) \cdot \Delta x.
$$
Hence, using the triangle inequality,
$$
\frac{\| f(x + \Delta x) - f(x) \|_Y}{\| \Delta x\|_X} = \frac{\| \nabla f (x) \cdot \Delta x + h(x + \Delta x) \cdot \Delta x \|_Y}{\| \Delta x\|_X} \leq \frac{\| \nabla f(x) \cdot \Delta x \|_Y}{\| \Delta x\|_X} + \| h(x + \Delta x) \|_{Y,X}.
$$
and
$$
\frac{\| f(x + \Delta x) - f(x) \|_Y}{\| \Delta x\|_X} = \frac{\| \nabla f (x) \cdot \Delta x + h(x + \Delta x) \cdot \Delta x \|_Y}{\| \Delta x\|_X} \geq \frac{\| \nabla f(x) \cdot \Delta x \|_Y}{\| \Delta x\|_X} - \| h(x + \Delta x) \|_{Y,X}.
$$
Furthermore, with $\xi$ taking the place of $\lambda \Delta x$ in the last step,
$$
\begin{align*}
\lim_{\delta\rightarrow 0} \sup_{0< \|\Delta x\|\leq \delta} \frac{\| \nabla f(x) \cdot \Delta x \|_Y}{\| \Delta x\|_X} \pm \| h(x + \Delta x) \|_{Y,X}
\end{align*} & = \lim_{\delta\rightarrow 0} \sup_{0< \|\Delta x\|\leq \delta} \frac{\| \nabla f(x) \cdot \Delta x \|_Y}{\| \Delta x\|_X}\\
& = \lim_{\delta\rightarrow 0} \sup_{0< \|\Delta x\|\leq \delta, 0 \neq \lambda \in \mathbb{R}} \frac{\| \nabla f(x) \cdot \lambda \Delta x \|_Y}{\| \lambda \Delta x\|_X}\\
& = \sup_{ \xi \neq 0} \frac{\| \nabla f(x) \cdot \xi \|_Y}{\| \xi \|_X} = \| \nabla f(x) \|_{Y,X}.
$$
We deduce
$$
\kappa_{abs}(x) = \lim_{\delta\rightarrow 0} \sup_{0<\|\Delta x\|\leq \delta} \frac{\| f(x + \Delta x) - f(x) \|_Y}{\| \Delta x\|_X} = \| \nabla f(x) \|_{Y,X}.
$$
Similarly,
$$
\kappa_{rel}(x) =
\frac{\| x \|_X}{\| f(x) \|_Y} \lim_{\delta\rightarrow 0} \sup_{0<\|\Delta x\|\leq \delta} \frac{\| f(x + \Delta x) - f(x) \|_Y}{\| \Delta x\|_X} = \frac{\| x \|_X}{\| f(x) \|_Y} \| \nabla f(x) \|_{Y,X}.
$$
:class: tip
Let $x \in \mathbb{R}^2$ and $f(x) = x_1 - x_2$. Compute the $\infty$-norm condition number $\kappa_{rel}(x)$ of $f(x)$. For what inputs is the condition number large?
Hint: Use the expression for the condition number of a differentiable function.
The Jacobian of $f$ is $J = \nabla f = \begin{pmatrix}1 & -1\end{pmatrix}$. Hence, $\|J\|_{\infty} = 2$ is the row-sum norm of $\begin{pmatrix}1 & -1\end{pmatrix}$. For the condition number, we obtain
$$
\kappa = \frac{\|x\|_\infty}{\|f(x)\|_\infty} \|J\|_\infty = \frac{2 \cdot \max\left\{|x_1|, |x_2|\right\}}{|x_1 - x_2|}.
$$
The condition number is large if $x_1 \approx x_2$. This reflects the issue of cancellation errors. Consider two numbers $x_1$ and $x_2$ that agree to the first 5 digits and each of them is accurate to 7 digits. The difference between the two numbers will only be accurate to 2 digits since the first 5 correct digits cancel each other out.
:class: tip
Let $f : D \subset X \to Y$ be differentiable and $\| x \|_X \neq 0$ and $\| f(x) \|_Y \neq 0$, $X$ and $Y$ being finite-dimensional real vector spaces. Show that
$$
\frac{\|\Delta f\|_Y}{\|f(x)\|_Y} \leq K_{rel}\cdot \frac{\|\Delta x\|_X}{\|x\|_X}
$$
and
$$
\frac{\|\Delta f\|_Y}{\|f(x)\|_Y} \leq \kappa_{rel}(x)\cdot \frac{\|\Delta x\|_X}{\|x\|_X} + \text{higher-order terms}
$$
We only show the second bound; the first follows from an inequality similar to $(*)$ below, with an adjusted set over which the supremum is taken. We multiply and divide by an auxiliary term: for $\|\Delta x\|_X > 0$,
$$
\frac{\|\Delta f\|_Y}{\|f(x)\|_Y} \leq \Bigl( \frac{\|\Delta f\|_Y}{\|f(x)\|_Y} \Big/ \frac{\|\Delta x\|_X}{\|x\|_X} \Bigr) \frac{\|\Delta x\|_X}{\|x\|_X} \stackrel{(*)}{\leq} \frac{\|\Delta x\|_X}{\|x\|_X} \sup_{0 < \| \Delta \tilde{x} \|_X \leq \| \Delta x \|_X} \Bigl( \frac{\| f(x + \Delta \tilde{x}) - f(x) \|_Y}{\| f(x) \|_Y} \Big/ \frac{\|\Delta \tilde{x}\|_X}{\|x\|_X} \Bigr)
$$
Because $\| f(x) \| \neq 0$ and $f$ is differentiable, we may use the arguments of the previous self-check question to show that the function
$$
g(\Delta x) := \sup_{0 < \| \Delta \tilde{x} \|_X \leq \| \Delta x \|_X} \Bigl( \frac{\| f(x + \Delta \tilde{x}) - f(x) \|_Y}{\| f(x) \|_Y} \Big/ \frac{\|\Delta \tilde{x}\|_X}{\|x\|_X} \Bigr)
$$
has the continuous extension
$$
g(0) = \frac{\| x \|_X}{\|f(x) \|_Y}\| \nabla f(x) \|_{Y,X} = \kappa_{rel}(x)
$$
at $\Delta x = 0$, meaning that there is a function $\tilde{h}: \mathbb{R}^n \to \mathbb{R}$ with $\lim_{\| \Delta x \| \to 0} \tilde{h}(\Delta x) = 0$ such that
$$
g(\Delta x) = \kappa_{rel}(x) + \tilde{h}(\Delta x).
$$
Now the result follows with
$$
\text{higher-order terms} = \tilde{h}(\Delta x) \frac{\|\Delta x\|_X}{\|x\|_X}.
$$