{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Numerical differentiation by finite diferences" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our goal is to explore the dependence of the numerical error on the size of perturbation of the independent variable in [finite difference (FD)](https://en.wikipedia.org/wiki/Numerical_differentiation#Finite_differences) numerical differentiation scheme." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As an example, consider a scalar function of a vector argument $f(\\mathbf x) = \\sum_{i=1}^n x_i^2$" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "f (generic function with 1 method)" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "f(x) = sum(xᵢ^2 for xᵢ in x) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's set some random vector $\\mathbf x\\in\\mathbb R^n$. First, we pick the dimension `n`" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "10" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "j = 1\n", "n = 10^j" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "10-element Array{Float64,1}:\n", " 0.8566789188220891\n", " 0.5816498901094684\n", " 0.6358555521874019\n", " 0.030924335530379432\n", " 0.9254341191018551\n", " 0.10640604323138736\n", " 0.9965536786537663\n", " 0.6642923827516156\n", " 0.1271063358389557\n", " 0.6479414279483935" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x₀ = rand(n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The function `f` evaluated at `x0` is" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4.215622236110449" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "f(x₀)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, in order to compute the gradient $\\nabla f$, we need to compute all the individual partial derivatives, the individual components of the vector. Let's now restrict ourselves just to one component, say, the first one, that is, let's compute $\\frac{\\partial f(\\mathbf x)}{\\partial x_1}$. In this simple example a formula can be written down upon inspection: $\\frac{\\partial f(\\mathbf x)}{\\partial x_1} = 2x_i$), but let's pretend that this answer is not available to us or we opt not to use it." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For reference, we compute the first component of the gradient using the \"exact\" formula above" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.7133578376441783" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dfdx_1_exact = 2*x₀[1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's now proceed to the numerical computation of this derivative." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Numerical computation using quotients (finite differences)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The perturbation $h$ of the first compoment of the vector $x$ is" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.001" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "k = 3\n", "h = 10.0^(-k)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Below we define a function for computing the first entry in the gradient (vector). Note that in defining the function we exploit the [multiple dispatch functionality of Julia](https://docs.julialang.org/en/v1/manual/methods/), thanks to which the function will compute with the floating point model of its input. That is, the input vector could be [IEEE-754 double-precision floating-point](https://en.wikipedia.org/wiki/Double-precision_floating-point_format) format, [IEEE-754 single-precision floating-point](https://en.wikipedia.org/wiki/Single-precision_floating-point_format) format (or perhaps even something else). " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "dfdx_1 (generic function with 1 method)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function dfdx_1(f,x₀::Vector{T},h::T) where T<:Real\n", " x = copy(x₀)\n", " x[1] = x[1]+h\n", " (f(x)-f(x₀))/h\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's compute the answer for the vector given (converted) into *IEEE single* format." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.7147063f0" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dfdx_1(f,Vector{Float32}(x₀),Float32(h))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's repeat the computation with the data in the *IEEE double* format (default for Julia)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.71435783764462" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dfdx_1(f,x₀,h)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Obviously, both answers differ from the correct one computed above " ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.7133578376441783" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dfdx_1_exact" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Dependence of the error on the length of the step" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's examine the error as a function of the size of the interval $h$" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "h_range = exp10.(range(-13, stop=-1, length=1000));" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "abs_er_64 = [abs((dfdx_1_exact - dfdx_1(f,x₀,h))) for h in h_range];\n", "abs_er_32 = [abs((dfdx_1_exact - dfdx_1(f,Vector{Float32}(x₀),Float32(h)))) for h in h_range];" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Plots.GRBackend()" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Plots\n", "gr()" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(h_range, abs_er_64,xaxis=:log, yaxis=:log, xlabel=\"h\", ylabel = \"|e|\", label = \"IEEE double\")\n", "plot!(h_range, abs_er_32,xaxis=:log, yaxis=:log, xlabel=\"h\", ylabel = \"|e|\", label = \"IEEE single\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we read the graph from right to left – as $h$ is getting smaller –, we observe that initially the error decreases for both 64- and 32-bit floating-point format, and it decreases at the same rate. The reason for this decrease is that the trunction error dominates here, and this error goes down with $h$.\n", "\n", "The major intended takeaway from this example is that this reduction of the error only takes place down to some $h$ below which no longer decreases; in fact, it actally increases as $h$ gets smaller. The reason is that for $h$ this small, the rounding errors dominate. Apparently, they start exhibitting themselves for smaller values with the 64-bit format than with 32-bit format. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is known from rigorous numerical analysis that the error in the case of the simple backward or forward finite difference approximation to a scalar derivative scales with $\\sqrt{\\epsilon}$, where $\\epsilon$ is the [machine epsilon](https://en.wikipedia.org/wiki/Machine_epsilon). Here we can expect even worse outcome as the dimension $n$ of the vector grows. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that $\\epsilon$ for double precision IEEE is" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.1102230246251565e-16" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "2^(-53)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "which is available in Julia through [eps()](https://docs.julialang.org/en/v1/manual/integers-and-floating-point-numbers/#Machine-epsilon) function with the type as the input argument (if no argument is given, it assumes `Float64`):" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.4901161193847656e-8" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sqrt(eps())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similarly, the 32-bit version is" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.00034526698f0" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sqrt(eps(Float32))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Apparently, these are roughly the \"cutt-off\" values of $h$." ] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "kernelspec": { "display_name": "Julia 1.5.3", "language": "julia", "name": "julia-1.5" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.5.3" } }, "nbformat": 4, "nbformat_minor": 2 }