{
"cells": [
{
"cell_type": "markdown",
"id": "4159be0e",
"metadata": {
"hideCode": false,
"hidePrompt": false,
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# 2022-01-19 Conditioning\n",
"\n",
"* Office hours: Monday 9-10pm, Tuesday 2-3pm, Thursday 2-3pm\n",
"* More on floating point\n",
"* Discuss Taylor Series activity\n",
"* Condition numbers\n",
"* Forward and backward error\n",
"* Computing volume of a polygon"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "eb781a7a",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [],
"source": [
"using Plots\n",
"default(linewidth=4)"
]
},
{
"cell_type": "markdown",
"id": "adc22a07",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Floating point representation is **relative**\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"id": "2022fd10",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"source": [
"Let $\\operatorname{fl}$ round to the nearest floating point number.\n",
"\n",
"$$ \\operatorname{fl}(x) = x (1 + \\epsilon), \\quad \\text{where} |\\epsilon| \\le \\epsilon_{\\text{machine}} $$\n",
"\n",
"This also means that the relative error in representing $x$ is small:\n",
"\n",
"$$ \\frac{|\\operatorname{fl}(x) - x|}{|x|} \\le \\epsilon_{\\text{machine}} $$"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "9007ab67",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(x -> (1 + x) - 1, xlims=(-1e-15, 1e-15))\n",
"plot!(x -> x)"
]
},
{
"cell_type": "markdown",
"id": "415173a8",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Exact arithmetic, correctly rounded"
]
},
{
"cell_type": "markdown",
"id": "3d1fcc95",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"source": [
"Take an elementary math operation $*$ (addition, subtraction, multiplication, division), and the discrete operation that our computers perform, $\\circledast$. Then\n",
"\n",
"$$x \\circledast y := \\operatorname{fl}(x * y)$$\n",
"\n",
"with a relative accuracy $\\epsilon_{\\text{machine}}$,\n",
"\n",
"$$ \\frac{|(x \\circledast y) - (x * y)|}{|x * y|} \\le \\epsilon_{\\text{machine}} . $$"
]
},
{
"cell_type": "markdown",
"id": "d5df8583",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"## Seems easy, how do operations compose?\n",
"\n",
"Is this true?\n",
"\n",
"$$ \\frac{\\Big\\lvert \\big((x \\circledast y) \\circledast z\\big) - \\big((x * y) * z\\big) \\Big\\rvert}{|(x * y) * z|} \\le^? \\epsilon_{\\text{machine}} $$"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "fb3c12ce",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f(x; y=1, z=-1) = (x+y)+z # The best arbitrary numbers are 0, 1, and -1\n",
"plot(x -> abs(f(x) - x)/abs(x), xlims=(-1e-15, 1e-15))"
]
},
{
"cell_type": "markdown",
"id": "9a0073c2",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Which operation caused the error?\n",
"\n",
"1. $\\texttt{tmp} = \\operatorname{fl}(x + 1)$\n",
"2. $\\operatorname{fl}(\\texttt{tmp} - 1)$\n",
"\n",
"Use Julia's [`BigFloat`](https://docs.julialang.org/en/v1/base/numbers/#BigFloats-and-BigInts)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "da0cfd0f",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(big(0.1)) = BigFloat\n",
"big(0.1) = 0.1000000000000000055511151231257827021181583404541015625\n",
"BigFloat(\".1\") = 0.1000000000000000000000000000000000000000000000000000000000000000000000000000002\n"
]
}
],
"source": [
"@show typeof(big(.1))\n",
"@show big(.1) # Or BigFloat(.1); parsed as Float64, then promoted\n",
"@show BigFloat(\".1\"); # Parse directly to BigFloat"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "22d3dadf",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.102230246251563524952071662733800140614440125894379682676737388538642032894338e-16"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tmp = 1e-15 + 1\n",
"tmp_big = big(1e-15) + 1 # Parse as Float64, then promote\n",
"abs(tmp - tmp_big) / abs(tmp_big)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "13055e31",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"r = tmp - 1\n",
"r_big = big(tmp) - 1\n",
"abs(r - r_big) / abs(r_big)"
]
},
{
"cell_type": "markdown",
"id": "2b931233",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Activity: [2022-01-12-taylor-series](https://classroom.github.com/a/VkPvGOgu)\n",
"\n",
"* Use Julia, Jupyter, Git\n",
"* Look at how fast series converge when taking only finitely many terms\n",
"* Explore instability, as is occuring for large negative `x` above, but not for the standard library `expm1`"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "84feb523",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [],
"source": [
"function myexp(x, k)\n",
" sum = 0\n",
" term = 1\n",
" n = 1\n",
" # modify so at most k terms are summed\n",
" while sum + term != sum\n",
" sum += term\n",
" term *= x / n\n",
" # YOUR CODE HERE\n",
" if n == k\n",
" break\n",
" end\n",
" n += 1\n",
" end\n",
" sum\n",
"end\n",
"\n",
"rel_error(x, k) = abs(myexp(x, k) - exp(x)) / exp(x)\n",
"ks = 2 .^ (0:10) # [1, 2, 4, ..., 1024];"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "fa9bbc4c",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(ks, k -> rel_error(-20, k), xscale=:log10, yscale=:log10)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "b85efacc",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(x -> rel_error(x, 1000) + 1e-17, xlims=(-20, 20), yscale=:log10)"
]
},
{
"cell_type": "markdown",
"id": "617909dc",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# What happened?\n",
"Let's look at the terms for positive and negative $x$"
]
},
{
"cell_type": "code",
"execution_count": 54,
"id": "41971b71",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"sum(expterms(x)) - exp(x) = 1.2092180894291739e-12\n",
"(sum(expterms(x)) - exp(x)) / exp(x) = 2.6634800885273225e-8\n"
]
},
{
"data": {
"text/plain": [
"51-element Vector{Float64}:\n",
" 1.0\n",
" -10.0\n",
" 50.0\n",
" -166.66666666666669\n",
" 416.66666666666674\n",
" -833.3333333333335\n",
" 1388.8888888888891\n",
" -1984.1269841269846\n",
" 2480.1587301587306\n",
" -2755.7319223985896\n",
" 2755.7319223985896\n",
" -2505.2108385441725\n",
" 2087.6756987868107\n",
" â‹®\n",
" -4.902469756513546e-8\n",
" 1.2256174391283866e-8\n",
" -2.9893108271424062e-9\n",
" 7.117406731291443e-10\n",
" -1.655210867742196e-10\n",
" 3.761842881232264e-11\n",
" -8.359650847182808e-12\n",
" 1.81731540156148e-12\n",
" -3.866628513960596e-13\n",
" 8.055476070751242e-14\n",
" -1.64397470831658e-14\n",
" 3.28794941663316e-15"
]
},
"execution_count": 54,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function expterms(x, k=50)\n",
" term = 1.\n",
" terms = [term]\n",
" for n in 1:k\n",
" term *= x / n\n",
" push!(terms, term)\n",
" end\n",
" terms\n",
"end\n",
"x = -10\n",
"@show sum(expterms(x)) - exp(x)\n",
"@show (sum(expterms(x)) - exp(x)) / exp(x)\n",
"expterms(x)"
]
},
{
"cell_type": "code",
"execution_count": 52,
"id": "5e1a6172",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"exp(-10) = 4.5399929762484854e-5\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 52,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@show exp(-10)\n",
"bar(expterms(-10))"
]
},
{
"cell_type": "markdown",
"id": "a2c8666b",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Conditioning\n",
"\n",
"> What sort of functions cause small errors to become big?\n",
"\n",
"Consider a function $f: X \\to Y$ and define the **absolute condition number**\n",
"$$ \\hat\\kappa = \\lim_{\\delta \\to 0} \\max_{|\\delta x| < \\delta} \\frac{|f(x + \\delta x) - f(x)|}{|\\delta x|} = \\max_{\\delta x} \\frac{|\\delta f|}{|\\delta x|}. $$\n",
"If $f$ is differentiable, then $\\hat\\kappa = |f'(x)|$.\n",
"\n",
"Floating point offers relative accuracy, so it's more useful to discuss **relative condition number**,\n",
"$$ \\kappa = \\max_{\\delta x} \\frac{|\\delta f|/|f|}{|\\delta x|/|x|}\n",
"= \\max_{\\delta x} \\Big[ \\frac{|\\delta f|/|\\delta x|}{|f| / |x|} \\Big] $$\n",
"or, if $f$ is differentiable,\n",
"$$ \\kappa = |f'(x)| \\frac{|x|}{|f|} . $$"
]
},
{
"cell_type": "markdown",
"id": "c33df973",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Condition numbers\n",
"\n",
"> $$ \\kappa = |f'(x)| \\frac{|x|}{|f|} $$"
]
},
{
"cell_type": "code",
"execution_count": 33,
"id": "e33212a9",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f(x) = x - 1; fp(x) = 1\n",
"plot(x -> abs(fp(x)) * abs(x) / abs(f(x)), xlims=(0, 2))"
]
},
{
"cell_type": "markdown",
"id": "cf70b0b6",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Back to $f(x) = e^x - 1$"
]
},
{
"cell_type": "code",
"execution_count": 55,
"id": "e3fbe77d",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 55,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f(x) = exp(x) - 1\n",
"fp(x) = exp(x)\n",
"plot(x -> abs(fp(x)) * abs(x) / abs(f(x)))"
]
},
{
"cell_type": "markdown",
"id": "b14994aa",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"## What does it mean?\n",
"\n",
"* The function $f(x) = e^x - 1$ is well-conditioned\n",
"* The function $f_1(x) = e^x$ is well-conditioned\n",
"* The function $f_2(x) = x - 1$ is ill-conditioned for $x \\approx 1$\n",
"\n",
"## The **algorithm** is unstable\n",
"\n",
"* `f(x) = exp(x) - 1` is unstable\n",
"* Algorithms are made from elementary operations\n",
"* Unstable algorithms do something ill-conditioned"
]
},
{
"cell_type": "markdown",
"id": "ba7e80f4",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# A stable algorithm"
]
},
{
"cell_type": "markdown",
"id": "ee9f8cc6",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"source": [
"We used the series expansion previously.\n",
"* accurate for small $x$\n",
"* less accurate for negative $x$ (see activity)\n",
"* we could use symmetry to fix\n",
"* inefficient because we have to sum lots of terms\n",
"\n",
"Standard math libraries define a more efficient stable variant, $\\texttt{expm1}(x) = e^x - 1$."
]
},
{
"cell_type": "code",
"execution_count": 35,
"id": "ccf0c2bb",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot([x -> exp(x) - 1,\n",
" x -> expm1(x)],\n",
" xlims = (-1e-15, 1e-15))"
]
},
{
"cell_type": "markdown",
"id": "5975cc87",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Another example $\\log(1 + x)$\n",
"\n",
"What is the condition number of $f(x) = \\log(1 + x)$ for $x \\approx 0$?"
]
},
{
"cell_type": "code",
"execution_count": 36,
"id": "8f2a60c3",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot([x -> log(1 + x),\n",
" x -> log1p(x)],\n",
" xlims = (-1e-15, 1e-15))"
]
},
{
"cell_type": "code",
"execution_count": 37,
"id": "4300fee4",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cond1(x) = abs(1/(1+x) * x / log1p(x))\n",
"cond2(x) = abs(1/x * x / log(x))\n",
"plot([cond1 cond2], xlims=(-1, 2), ylims=(0, 100))"
]
},
{
"cell_type": "markdown",
"id": "ca3502f6",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Reliable = well-conditioned and stable\n",
"\n",
"## Mathematical functions $f(x)$ can be ill-conditioned (big $\\kappa$)\n",
"* Modeling is how we turn an abstract question into a mathematical function\n",
"* We want well-conditioned models (small $\\kappa$)\n",
"* Some systems are intrinsically sensitive: fracture, chaotic systems, combustion\n",
"\n",
"## Algorithms `f(x)` can be unstable\n",
"* Unreliable, though sometimes practical\n",
"* Unstable algorithms are constructed from ill-conditioned parts"
]
},
{
"cell_type": "markdown",
"id": "07cd924a",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# An ill-conditioned problem from Paul Olum\n",
"\n",
"From [Surely You're Joking, Mr. Feynman](https://sistemas.fciencias.unam.mx/%7Ecompcuantica/RICHARD%20P.%20FEYNMAN-SURELY%20YOU%27RE%20JOKING%20MR.%20FEYNMAN.PDF) (page 113)\n",
"\n",
"> So Paul is walking past the lunch place and these guys are all excited. \"Hey, \n",
"Paul!\" they call out. \"Feynman's terrific! We give him a problem that can be stated in ten \n",
"seconds, and in a minute he gets the answer to 10 percent. Why don't you give him one?\" \n",
"Without hardly stopping, he says, \"The tangent of 10 to the 100th.\" \n",
"I was sunk: you have to divide by pi to 100 decimal places! It was hopeless."
]
},
{
"cell_type": "markdown",
"id": "6ecb7b9c",
"metadata": {
"cell_style": "split"
},
"source": [
"What's the condition number?\n",
"\n",
"$$ \\kappa = |f'(x)| \\frac{|x|}{|f|} $$\n",
"\n",
"* $f(x) = \\tan x$\n",
"* $f'(x) = 1 + \\tan^2 x$\n",
"\n",
"$$ \\kappa = \\lvert x \\rvert \\Bigl( \\lvert \\tan x \\rvert + \\bigl\\lvert \\frac{1}{\\tan x} \\bigr\\rvert \\Bigr)$$"
]
},
{
"cell_type": "code",
"execution_count": 56,
"id": "3c52979f",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"-0.4116229628832498"
]
},
"execution_count": 56,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tan(1e100)"
]
},
{
"cell_type": "code",
"execution_count": 63,
"id": "f20ecd97",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"0.4012319619908143541857543436532949583238702611292440683194415381168718098221194"
]
},
"execution_count": 63,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tan(BigFloat(\"1e100\", precision=500))"
]
},
{
"cell_type": "markdown",
"id": "9a15e2e3",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Go find some functions...\n",
"\n",
"* Find a function $f(x)$ that models something you're interested in\n",
"* Plot its condition number $\\kappa$ as a function of $x$\n",
"* Plot the relative error (using single or double precision; compare using Julia's `big`)\n",
"* Is the relative error ever much bigger than $\\kappa \\epsilon_{\\text{machine}}$?\n",
"* Can you find what caused the instability?\n",
"* Share on Zulip\n",
"\n",
"## Further reading: [FNC Introduction](https://fncbook.github.io/fnc/intro/overview.html)"
]
}
],
"metadata": {
"celltoolbar": "Slideshow",
"hide_code_all_hidden": false,
"kernelspec": {
"display_name": "Julia 1.7.1",
"language": "julia",
"name": "julia-1.7"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.7.1"
},
"rise": {
"enable_chalkboard": true
}
},
"nbformat": 4,
"nbformat_minor": 5
}