Finite Differences

Resources

Numerical Analysis Notes
https://tobydriscoll.net/fnc-julia/localapprox/fd-converge.html
Burden, Richard L., and J. Douglas Faires. 2010. Numerical Analysis. 9th ed. Cengage Learning.

Main Idea

TODO

Error Analysis

For a finite difference method approximating a $d$ th order derivative, with $k$ th order convergence, the error depends on $n_{x}$ as

\underset{Truncation Error}{\underset{⏟}{C_{1} n_{x}^{- k}}} + \underset{Floating Point Error}{\underset{⏟}{C_{2} n_{x}^{d}}} .

Or equivalently,

\underset{Truncation Error}{\underset{⏟}{C_{1} Δ x^{k}}} + \underset{Floating Point Error}{\underset{⏟}{C_{2} Δ x^{- d}}} .

$C_{1}$ depends on a bound of the $(d + k)$ th derivative, while $C_{2}$ depends linearly on the machine precision, and the conditioning of the function. Setting the two errors equal to one another gives the lowest error, revealing how the best mesh scales:

\begin{aligned} C_{1} n_{x}^{- k} & = c_{2} ε_{mach} \cdot n_{x}^{d} \\ \frac{C_{1}}{c_{2}} \cdot ε_{mach}^{- 1} & = n_{x}^{d + k} \\ n_{x} & \sim ε_{mach}^{- 1 / (d + k)} \end{aligned}

E.g., for a first-order derivative, the optimal mesh is

n_{x} \approx ε_{mach}^{- 1 / (k + 1)}, or Δ x = ε_{mach}^{1 / (k + 1)},

and in the more general case,

n_{x} \approx ε_{mach}^{- 1 / (k + d)}, or Δ x \approx ε_{mach}^{1 / (k + d)},

Changing from 32 bit to 64 bit Floating Point representations is like multiplying by this factor (the constants should cancel out). See the table below for a break-down of the change to the optimal $n_{x}$ that comes from this switch. Just squaring $n_{x}$ would unjustly include the constants twice, but still shows the general trend.

Floating Point Error Details

We wish to estimate $f^{'} (x)$ . Due to floating point error, we assume

fl (x) = (1 + ε) x,

where $| ε | \leq ε_{machine precision} \approx 10^{- 16} .$ Denoting the true finite difference (no floating point error) as

δ (x) = \frac{f (x + h) - f (x)}{h} .

Thus,

f^{'} (x) = δ (x) + \underset{O (h)}{\underset{⏟}{τ_{f} (h)}} .

Recalling the relative Condition Number,

κ_{f} (x) = | \frac{x f^{'} (x)}{f (x)} |,

and using a finite difference to approximate the numerator,

κ_{f} (x) \cdot | x | \approx | \frac{f (x + h) - f (x)}{h} | .

The resulting error is a combination of truncation and floating point error:

e (h) = | \partial_{x} f (x) - Δ_{x}^{h} f (x) | \leq \frac{ε}{h} + τ_{f} (h),

where for instance $τ_{h} = O (h^{2})$ for the 3-point centered method.

Taking $τ_{h} = (M h^{2}) / 6$ , where $M$ bounds $f^{‴}$ , the minimum error occurs at

h = {(\frac{3 ε}{M})}^{1 / 3} .

Practical Floating Point Error

Float32 vs. Float64 / double

How much does the best case error increase by switching from Float32 to Float64? Deep Learning commonly relies on Float32 (or even Automatic Mixed Precision) precision to reduce memory requirements or allow more parameters. This decision often trickles into scientific machine learning, but at what accuracy cost? For a given precision, what is the optimal $n_{x}$ ?

We can continue, and rearrange $h = (c ε)^{1 / 3}$ , where $c = κ_{f} (x) f (x)$ . Taking $c = 1$ , we have

ε = 10^{- 16} \to h = 10^{- 16 / 3} = 10^{- 5. \bar{3}} \approx 4.64158883 \cdot 10^{- 6}

ε = 10^{- 8} \to h = 10^{- 8 / 3} = 10^{- 2. \bar{6}} \approx 2.15443469 \cdot 10^{- 3}

For $[- 1, 1]$ , this would mean $n_{x} = 430886$ vs. $n_{x} = 928$ .

n_{x} = \frac{2}{h} = 2 \cdot c^{- 1 / 3} ε^{- 1 / 3}

Let's suppose with $ε_{f} = 10^{- 8}$ , $n_{x}^{f} = 100$ . Then, what would $c$ be, and what would $n_{x}^{d}$ be for $ε_{d} = 10^{- 16}$ ?

\begin{aligned} h^{f} & = 2 / n_{x}^{f} \\ ε_{d} & = ε_{f}^{2} \\ c & = {(\frac{2}{n_{x}^{f}})}^{3} \cdot \frac{1}{ε_{f}} = {(\frac{2}{n_{x}^{d}})}^{3} \cdot \frac{1}{ε_{d}} \\ \frac{8}{(n_{x}^{f})^{3}} \cdot \frac{1}{ε_{f}} & = \frac{8}{(n_{x}^{d})^{3}} \cdot \frac{1}{ε_{f}^{2}} \\ \frac{1}{(n_{x}^{f})^{3}} & = \frac{1}{(n_{x}^{d})^{3}} \cdot \frac{1}{ε_{f}} \\ (n_{x}^{f})^{3} & = (n_{x}^{d})^{3} \cdot ε_{f} \\ n_{x}^{f} & = n_{x}^{d} \cdot (ε_{f})^{1 / 3} \approx n_{x}^{d} \cdot \frac{1}{500} \\ n_{x}^{d} & = n_{x}^{f} \cdot (ε_{f})^{- 1 / 3} \approx n_{x}^{f} \cdot 500 \end{aligned}

In other words, by using double precision, we can use $\times 500$ as many points in our mesh as using single precision. On the other hand, if $ε_{f} = 10^{- 5.3}$ , this factor is $\times 58.4$ .

Example results for $\sin (x)$ on $[- 1, 1]$ :

Derivative Order	Truncation Order	Predicted Factor	Observed Factor	32bit $n_{x}$ , $n_{x}^{f}$
1	2	500	924	727
2	2	100	150	126
3	2	40	42	92
1	4	40	41	92
2	4	22	29	41

Here, these "Factors" are $c_{fac}$ such that $n_{x}^{d} = c_{fac} \cdot n_{x}^{f}$ , which we can either predict analytically based on the error form and the values of $ε$ ("Predicted"). Alternatively we can observe these factors by testing a bunch of $n_{x}$ on a smooth example problem and compare the optimal $n_{x}$ between the different precisions ("Observed").