Math Fundamentals

Numerical Methods

David Mayerich

Scalable Tissue Imaging and Modeling (STIM) Laboratory

Department of Electrical and Computer Engineering

Cullen College of Engineering

University of Houston

David Mayerich

STIM Laboratory, University of Houston

Taylor Series

Maclaurin Series

Error Terms

Approximating Complex Functions

David Mayerich

STIM Laboratory, University of Houston

Taylor Series

A special function has a standardized name and notation:
\(\quad \sin(x) \quad\quad e^x \quad\quad \tan^{-1}(x) \quad\quad \ln(x) \quad\quad \log_2(x) \quad\quad \text{erf}(x) \quad\quad J_n^{(1)}(x) \)
Digital processors perform basic numerical operations (ex. addition, multiplication)
Special functions require specific algorithms to evaluate
Start with one you probably know: Taylor expansions

David Mayerich

STIM Laboratory, University of Houston

\ln(1+x) = x - \frac{x^2}{2} + \frac{x^3}{3} - \frac{x^4}{3} \cdots

\sin(x) = x - \frac{x^3}{6} + \frac{x^5}{120} \cdots

But now we have infinite summations...

\sin(x) = \sum_{n=1}^\infty \frac{(-1)^n x^{1+2n}}{(1+2n)!}

\ln(1+x) = -\sum_{n=1}^\infty \frac{(-1)^n x^{n}}{n}

Consists of basic mathematical operations executed in a loop

Maclaurin Series

Start with expansions about zero
Creating an n-order Maclaurin polynomial requires differentiating a function n times

David Mayerich

STIM Laboratory, University of Houston

m(x)=f(0) + \frac{f^{'}(0)}{1!}x + \frac{f^{''}(0)}{2!}x^2 + \frac{f^{(3)}}{3!}x^3 + \frac{f^{(4)}(0)}{4!}x^4 \cdots

m(x) = \sum^{\infty}_{n=0} \frac{f^{(n)}(0)}{n!}x^n

How does this work?

m(0) = f(0) + \frac{f^{'}(0)}{1}\cdot 0 + \frac{f^{''}(0)}{2} \cdot 0 \cdots

m^{'}(x) = 1 \cdot \frac{f^{'}(0)}{1}\cdot x^{1-1} + 2 \cdot \frac{f^{''}(0)}{2} \cdot x^{2-1} \cdots

m^{'}(0) = f^{'}(0) + 2 \cdot \frac{f^{''}(0)}{2} \cdot 0 \cdots

Expansion of a Polynomial

David Mayerich

STIM Laboratory, University of Houston

p(x)=4x^3 + 2x^2 + x + 3

m(x)=p(0) + \frac{p'(0)}{1!}x + \frac{p''(0)}{2!}x^2 + \frac{p^{(3)}(0)}{3!}x^3 + \frac{p^{(4)}(0)}{4!}x^4 + \cdots

m(x)=3 + \frac{1}{1!}x + \frac{4}{2!}x^2 + \frac{24}{3!}x^3 + \frac{0}{4!}x^4 + \cdots

m(x)=3 + x + 2x^2 + 4x^3 + 0x^4 + \cdots

p(0)=3

p'(0)= 12(0)^2 + 4(0) + 1 = 1

p''(0)=24(0) + 4 = 4

p^{(3)}(0)=24

p^{(4)}(0)=0

m(x)=p(x)

Truncation Error

The Taylor expansion for a degree d polynomial has at most d non-zero terms:

David Mayerich

STIM Laboratory, University of Houston

Many functions have infinite approximations:

e^x = \sum_{n=0}^{\infty} \frac{x^n}{n!} = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \cdots

\cos(x) = \sum_{n=0}^{\infty} \frac{(-1)^n}{(2n)!} x^{2n} = 1 - \frac{x^2}{2!} + \frac{x^4}{4!} - \cdots

Calculating these functions requires truncating to a finite number of terms:

\sqrt{1+x} = 1 + \frac{1}{2}x - \frac{1}{8}x^2 + \frac{1}{16}x^3 + \epsilon

p(x)=4x^3 + 2x^2 + x + 3

m(x)=3 + x + 2x^2 + 4x^3 + 0x^4 + \cdots

we want this error to be small

\sqrt{1+x} \approx 1 + \frac{1}{2}x - \frac{1}{8}x^2 + \frac{1}{16}x^3

Remainder Term

The error term can be expressed analytically using the Lagrange remainder:

David Mayerich

STIM Laboratory, University of Houston

f(x)=f(0) + \frac{f^{'}(0)}{1!}x + \frac{f^{''}(0)}{2!}x^2 + \frac{f^{(3)}}{3!}x^3 + \frac{f^{(4)}(0)}{4!}x^4 \cdots

f(x) \approx f(0) + \frac{f^{'}(0)}{1!}x + \frac{f^{''}(0)}{2!}x^2 + \cdots + \frac{f^{(k)}(0)}{k!}x^k = \sum_{n=0}^{k} \frac{f^{(n)}(0)}{n!}(x)^n

f(x) = \sum_{n=0}^{k} \frac{f^{(n)}(0)}{n!}(x)^n + \frac{f^{(k+1)}(z)}{(k+1)!}x^{k+1}

We approximate the function by taking k terms in the series:

where z is some value \(0\leq z \leq x\)

What can we know about this error?

Lagrange Remainder Term

How will this term behave as the number of terms increases?
- Consider: factorial vs. exponential growth, numerator as derivative increases
How does this term change with \(x\) and \(k\)?
Consider the error with the following functions:

David Mayerich

STIM Laboratory, University of Houston

\epsilon = \frac{f^{(k+1)}(z)}{(k+1)!}x^{k+1}

f(x) = \sin(x)

f(x) = x^5 + 1

f(x) = e^x

Taylor Series

A Taylor series shifts the expansion away from zero

David Mayerich

STIM Laboratory, University of Houston

m(x)=f(0) + \frac{f'(0)}{1!}x + \frac{f''(0)}{2!}x^2 + \frac{f^{(3)}(0)}{3!}x^3 + \frac{f^{(4)}(0)}{4!}x^4 + \cdots

t(x)=f(x_0) + \frac{f'(x_0)}{1!}(x-x_0) + \frac{f''(x_0)}{2!}(x-x_0)^2 + \frac{f^{(3)}(x_0)}{3!}(x-x_0)^3 + \cdots

x_0=0:

How can this benefit us?

x_0=-1

x_0=0

x_0=1

Sets and Spaces

Sets of Numbers

Cartesian Products

Tensors

David Mayerich

STIM Laboratory, University of Houston

Number Sets

Natural Numbers (can start from 0 or 1)
Integers (includes negative numbers)
Rational Numbers
- fractions of integers
Real Numbers
Complex Numbers

David Mayerich

STIM Laboratory, University of Houston

\mathbb{N} \in \left\{ 0, 1, 2, \cdots \right\}

\mathbb{Z} \in \left\{ \cdots -2, -1 0, 1, 2, \cdots \right\}

\mathbb{Q}

\mathbb{C} \in a + bi \quad \text{where} \quad a,b \in \mathbb{R} \quad \text{and} \quad i^2 = -1

\mathbb{I}

\mathbb{C}

\mathbb{N}

\mathbb{Z}

\mathbb{Q}

\mathbb{R}

-1

-\frac{1}{2}

-\sqrt{2}

-\pi

\pi

\frac{1}{2}

\sqrt{2}

\mathbb{R} \in

Building Sets

Set Builder Notation
- variable: potential member of a set
- predicate: condition for membership
Common conditional operations:
Real numbers in the range
A domain is commonly specified with the variable:
Evaluate the following sets:

David Mayerich

STIM Laboratory, University of Houston

\left\{ x | x > 0 \right\}

\left\{x | x \in \mathbb{R} \land x \geq 0 \land x < 5 \right\}

a \land b

logical and

a \lor b

logical or

\exist y

"there exists"

[0, 5)

\left\{x \in \mathbb{R} | x \geq 0 \land x < 5 \right\}

\left\{x \in \mathbb{R} | x=x^2 \right\}

\left\{n \in \mathbb{Z} | (\exists k \in \mathbb{Z})[n=2k] \right\}

\left\{x \in \mathbb{R} | (\exists p \in \mathbb{Z})(\exists q \in \mathbb{Z})\left[x\neq 0 \land x = \frac{p}{q} \right] \right\}

\rightarrow \{0,1\}

\rightarrow \text{even integers}

\rightarrow \mathbb{Q}

variable

predicate

Cartesian Product and n-Tuples

David Mayerich

STIM Laboratory, University of Houston

\text{for two sets } A \text{ and } B \text{:} \quad A \times B = \left\{ (a, b) | a\in A \land b \in B \right\}

\mathbb{R} \times \mathbb{R} = \left\{ (x, y) | x\in \mathbb{R} \land y \in \mathbb{R} \right\} = \mathbb{R}^2

\mathbb{C} \times \mathbb{C} \times \mathbb{C} = \left\{ (x, y, z) | x\in \mathbb{C} \land y \in \mathbb{C} \land z \in \mathbb{C} \right\} = \mathbb{C}^3

\begin{bmatrix} x \in \mathbb{R}\\ y \in \mathbb{R} \end{bmatrix}

\begin{bmatrix} x \in \mathbb{C}\\ y \in \mathbb{C}\\ z \in \mathbb{C} \end{bmatrix}

\mathbb{R}^2 \times \mathbb{R}^2 = \left\{ (x, y) \in \mathbb{R}^2 \right\} \times \left\{ (u, v) \in \mathbb{R}^2 \right\}

\begin{bmatrix} x \in \mathbb{R} & u \in \mathbb{R}\\ y \in \mathbb{R} & v \in \mathbb{R} \end{bmatrix}

\mathbb{R}^3 \times \mathbb{R}^3 = \left\{ (x, y, z) \in \mathbb{R} \right\} \times \left\{ (u, v, w) \in \mathbb{R} \right\} = \left(\mathbb{R}^3\right)^2 = \mathbb{R}^{3 \times 2}

\begin{bmatrix} x \in \mathbb{R} & u \in \mathbb{R}\\ y \in \mathbb{R} & v \in \mathbb{R}\\ z \in \mathbb{R} & w \in \mathbb{R} \end{bmatrix}

2-tuple

3-tuple

= \left\{ (x, y), (u,v) | (x, y) \in \mathbb{R}^2 \land (u, v) \in \mathbb{R}^2 \right\}

= \left(\mathbb{R}^2\right)^2=\mathbb{R}^{2 \times 2}

Tensor and Outer Products

The tensor product is a Cartesian product

that also satisfies the following distributive properties:
For sets of numbers these conditions lead to the "outer product" in linear algebra:

David Mayerich

STIM Laboratory, University of Houston

A \otimes B = A \times B

(a + c) \otimes b = a \otimes b + c \otimes b

a \otimes (b + c) = a \otimes b + b \otimes c

\lambda (a \otimes b) = (\lambda a) \otimes b = a \otimes (\lambda b)

\begin{bmatrix} x \\ y \end{bmatrix} \in \mathbb{R}^2

\begin{bmatrix} u \\ v \end{bmatrix} \in \mathbb{R}^2

\begin{bmatrix} x \\ y \end{bmatrix} \otimes \begin{bmatrix} u \\ v \end{bmatrix}

= \begin{bmatrix} x \\ y \end{bmatrix} \begin{bmatrix} u & v \end{bmatrix}

= \begin{bmatrix} xu & xv \\ yu & yv \end{bmatrix} \in \mathbb{R}^{2 \times 2}

Tensors and Arrays

A set produced using tensor products is a tensor
The number of tensor products defines the order of the tensor
0-order tensors are scalar values:
1-order tensors are vectors
2-order tensors are matrices
Higher-order (3+) tensors exist
- Mathematics is just extended
- Machine learning, image processing, physics

David Mayerich

STIM Laboratory, University of Houston

a\in \mathbb{R} \quad b\in \mathbb{C}

\mathbf{v} = \begin{bmatrix} v_1\\ v_2\\ \vdots\\ v_n \end{bmatrix} \in \mathbb{R}^n

\mathbf{A} = \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1m}\\ a_{21} & a_{22} & \cdots & a_{2m}\\ \vdots & \vdots & \ddots & \vdots\\ a_{n1} & a_{n2} & \cdots & a_{nm}\\ \end{bmatrix} \in \mathbb{R}^{n\times m}

Manipulating Tensors

Representing tensors in algorithms:

David Mayerich

STIM Laboratory, University of Houston

Use notation to constrain your requirements when designing algorithms:

what is the vector space of B and D?

(\mathbf{A} + \mathbf{B}) \mathbf{C} = \mathbf{D} \quad \text{where} \quad \mathbf{A} \in \mathbb{R}^{3\times 5} \quad \text{and} \quad \mathbf{C} \in \mathbb{C}^{5\times 4}

\mathbf{B} \in ?^{3 \times 5}

\mathbf{D} \in \mathbb{C}^{3\times 4}

\left( \mathbb{R}^{3\times 5} + ? \right) \mathbb{C}^{5\times 4} = ?

\left( \mathbb{R}^{3\times 5} + ?^{3\times 5} \right) \mathbb{C}^{5\times 4} = \mathbb{C}^{3\times 4}

We don't know the set of B, but we do know it's size.

We know the sum is at least

\mathbb{C}^{3\times 5}

\mathbf{a} \in \mathbb{R}^{5}

float a[5];
float* a = (float*) malloc(5 * sizeof(float));
float* a = new float[5];

\mathbf{B} \in \mathbb{C}^{3\times 4}

float B[3][4][2];
float* B = (float*) malloc(3 * 4 * 2 * sizeof(float));
std::complex<float>* = new std::complex<float>[3 * 4];

a = numpy.zeros((5, 1))

a = numpy.zeros((3, 4), dtype=numpy.complex64)

C/C++

Python

Functions

Continuous Functions

Domain and Range Sets

Discrete Functions

David Mayerich

STIM Laboratory, University of Houston

Continuous and Discrete

Functions provide a mapping from one space to another
Consider the Bessel function (often used in EM scattering):
This function is difficult to calculate for any \(x\) and \(n\)
- Evaluate the series above and truncate
- Use another approximation:
Pre-compute a discrete array and access the values as needed: \(J[n_i, x_j]\)
- I will use (...) to represent to continuous functions and [...] for discrete functions
- The values \( n_i \in \mathbb{N} \) and \(x_j \in \mathbb{N}\) are indices calculated from \(n\) and \(x\)

David Mayerich

STIM Laboratory, University of Houston

J(n, x) = \sum_{m=0}^{\infty}\frac{(-1)^m\left(\frac{x}{2} \right)^{n+2m}}{m!(n+m)!}

\text{where} \quad J\in\mathbb{R} \quad n \in \mathbb{N} \quad x \in \mathbb{R}

J(n, x) \approx \sqrt{\frac{2}{\pi x}} \cos\left( x - \left[\frac{n}{2} + \frac{1}{4} \right] \right)

Discrete Function Examples

How can we define the set for a 1024 x 768 color image?
The "theoretical" image can be expressed as a continuous function:
\( \mathbf{m}(\mathbf{x})\in \mathbb{R}^3\) where \(\mathbf{x}\in\mathbb{R}^2\)
- Each point returns a 3-tuple (representing red, green, and blue)
- This is equivalent to: \(\mathbf{m}(x, y)\in\mathbb{R}^3\) where \( x, y \in \mathbb{R}\)
In computer memory the image will be represented as a series of samples:
\( \mathbf{m}[\mathbf{u}]\in \mathbb{R}^3\) where \(\mathbf{u}\in\mathbb{N}^2\)

David Mayerich

STIM Laboratory, University of Houston

Hilbert Spaces

Complete Sets

Inner Products

David Mayerich

STIM Laboratory, University of Houston

Complete Sets

A set \(\mathbb{M}\) is incomplete if a series of elements in \(\mathbb{M}\) converges to a result outside of \(\mathbb{M}\)

David Mayerich

STIM Laboratory, University of Houston

\frac{\pi}{4} = 1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} \cdots = \sum_{n=0}^{\infty} \frac{(-1)^n}{2n + 1}

\tan^{-1}(x) = x - \frac{x^3}{3} + \frac{x^5}{5} - \frac{x^7}{7} \cdots = \sum_{n=0}^{\infty} \frac{(-1)^n x^{n+1}}{2n + 1}

or calculating the arc-tangent:

A complete space includes all points in a series and the values they converge to

(\(\mathbb{R}\) and \(\mathbb{C}\) are complete, \(\mathbb{Q}\) is not)

Consider the Lebniz formula for \(\pi\):

\mathbb{Q}

\mathbb{R}

(if \(x \in \mathbb{Z}\))

Inner Products

An inner product \(\langle\mathbf{a} , \mathbf{b}\rangle\) is an operation across some set \(\mathbb{M}\) with the following properties:

maps two vector values to a scalar: \(\langle\mathbf{a},\mathbf{b}\rangle=c\) where \(\mathbf{a}\in\mathbb{M}^n\), \(\mathbf{b}\in\mathbb{M}^n\), and \(c\in\mathbb{M}\)
conjugate symmetric: \(\langle\mathbf{a},\mathbf{b}\rangle = \overline{\langle\mathbf{b},\mathbf{a}\rangle}\), where \(\overline{x}\) is the conjugate transpose of \(x\)
positive definite: \(\langle\mathbf{a},\mathbf{a}\rangle >0\), unless \(\mathbf{a} = \mathbf{0}\) in which case \(\langle \mathbf{a},\mathbf{a}\rangle = 0\)
linear in the first argument: \(\langle \alpha\mathbf{a} +\alpha\mathbf{b}, \mathbf{c}\rangle = \alpha\langle \mathbf{a}, \mathbf{c}\rangle + \alpha \langle \mathbf{b}, \mathbf{c}\rangle\)

David Mayerich

STIM Laboratory, University of Houston

The inner product in \(\mathbb{R}\) and \(\mathbb{C}\) is usually defined as the "dot product":

\langle \mathbf{x} , \mathbf{y} \rangle = \mathbf{x} \cdot \mathbf{y} = \mathbf{x}^T \mathbf{y} = \sum_{i=1}^{n} x_i y_i \quad \text{when} \quad \mathbf{x},\mathbf{y}\in\mathbb{R}^n

\langle \mathbf{x} , \mathbf{y} \rangle = \mathbf{y}^\dag \mathbf{x} = \sum_{i=1}^{n} x_i \overline{y}_i \quad \text{when} \quad \mathbf{x},\mathbf{y}\in\mathbb{C}^n

Hilbert Spaces

A Hilbert space is a set that is:

complete
has an inner product

David Mayerich

STIM Laboratory, University of Houston

Real and complex vector spaces \(\mathbb{R}^n\) and \(\mathbb{C}^n\) are the most common examples.

Another example includes the set of all functions (ex. \(f(x)\in\mathbb{C}\)). The inner product between two functions \(f(x)\) and \(g(x)\) is:

\langle f , g \rangle = \int_{a}^{b} f(t) \overline{g}(t) dt

which is a continuous extension of the vector dot product:

\langle \mathbf{x} , \mathbf{y} \rangle = \mathbf{y}^\dag \mathbf{x} = \sum_{i=1}^{n} x_i \overline{y}_i

Discussion

Does the set of all polynomials \(\mathbb{P} = \{x^0, x^1, x^2, \cdots \}\) form a Hilbert space?
What is the upper bound (worst-case) error if we approximate \(e^x\) using a 5th-order Maclaurin series in the range \(x\in [0, 1]\)?
Assume that we are running a simulation to calculate the electric field \(\mathbf{E}(\mathbf{r})\), which describes the scattered field at a position \(\mathbf{r}\).
- What is the Hilbert space for \(\mathbf{E}\)?
- What is the Hilbert space for \(\mathbf{r}\)?
- How would these change when running the simulation on a discrete grid?

David Mayerich

STIM Laboratory, University of Houston

\(\mathbf{r}\in \mathbb{R}^3\)

\(\mathbf{E}(\mathbf{r},t)\in \mathbb{R}^3\)

\(\mathbf{E}(\mathbf{r})\in \mathbb{C}^3\)

\epsilon = \frac{f^{(k+1)}(z)}{(k+1)!}x^{k+1}

Lagrange remainder: