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}
0
-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:
A.2 Math Fundamentals
By STIM Laboratory
