Numerical Differentiation in Coding: The Pythonic Way
Have you had problems coding the differential value of a function f(x)
? Do you need a functional approach that can automate differentiation for you? If the answer to either of these queries is a yes, then this blog post is definitely meant for you.
I have provided the Google Colaboratory Notebook and the Github Repository Num_Diff so you may follow along using these utilities for a better understanding of the topic.
Why Numerical Differentiation?
Consider the following problem :
Before writing some formal code, import the math module
in Python:
import math
Now, you might construct the objective function as follows:
def f(x):
return math.asin(math.sqrt(1-math.pow(x,2)))
For the differential function f’(x)
, a Typical Mathematician would thereafter go on to manually differentiate the function using lengthy and complex chain rules, that might look as follows and then code it:
def differentialfx(x):
return (-x)/(abs(x) * math.sqrt(1-math.pow(x,2)))
Scary, eh? Well, we’re not here for that, we will try to automate the differentiation operator(d/dx
) as a whole.
Concept of Limits, Derivatives, and Approximations
The limit of a function at a point
a
in its domain (if it exists) is the value that the function approaches as its argument approachesa
.
We will use the above definition to obtain the Left & Right Hand Derivatives, respectively as follows:
For better approximation often the average of the two is considered the best estimate of the differential at x=a
. Mathematically, what we are trying to compute is:
Hence to get the differential value of f(x) at x=a [f’(x=a)
] we will use the above formula, substituting a very small positive number {say 10^(-6)
} for h
. Let’s dive right into coding it :
def differentiate(func,a,h=1e-6):
"""
Returns a derivative of the passed function f(x) at a given value
a using the concept of Right Hand and Left Hand derivative
approximation.
f = univariate function f(x)
a = value at which derivative should be calculated f'(x=a)
h = the tolerance h, which is a value very close to zero but not
zero.
"""
return (func(a+h)-func(a-h))/(2*h)
Results with Verification
Now let us find out some differentials using the above-defined scripts for 0<x<1
. The approach we will use here is simple :
- Loop through
num
from0.01 to 0.99
- Calculate differential using hard-coded method & numerical differentiation
- Know if the value of the two computed values are close to each other or not.
- If the values aren’t close, store and display the value.
- Inspect the individual values where the results were not close.
DIFFS = []
for i in range(1,100):
num = i/100
approx_estimated_value = differentiate(f,num)
actual_true_value = differentialfx(num)
val = math.isclose(approx_estimated_value,
actual_true_value,rel_tol=1e-9)
if val == False:
print("Suitable Difference Found at",num)
DIFFS.append(num)
Output Obtained :
Suitable Difference Found at 0.02
Suitable Difference Found at 0.99
Note: These differences are occurring because of approximation errors. If we wish to check till
n
digits after the decimal, these differences can then be avoided by usingn=6
(or a lower value) inmath.isclose(approx_estimated_value, actual_true_value,rel_tol=1e-n)
. In this case no drastic differences are found.
Inspect the Differences :
for each_num in DIFFS:
print("Values at",each_num,":",differentiate(f,each_num,h=1e-6),
",",differentialfx(each_num))
Output:
Values at 0.02 : -1.0002000616626816 , -1.000200060020007
Values at 0.99 : -7.088812059172223 , -7.088812050083354
The values only differ after 8 digits or so, hence running the script with math.isclose(approx_estimated_value, actual_true_value,rel_tol=1e-8)
shows that there are no differences.
Higher Level Information
PyTorch (a Python deep learning module) has AutoGrad Features that employ Chain Rule Mechanisms of Differential Calculus using complex tree-like structures (graphs) that perform the same in a more efficient and faster way:
import torch
x = torch.autograd.Variable(torch.Tensor([0.5]),requires_grad=True)
def fnew(x):
return torch.asin(torch.sqrt(1-torch.pow(x,2)))
y = fnew(x)
y.backward()
print(float(x.grad))
However, the implementation from the math
module will not work whilst working with PyTorch. We have to use the torch
module to use autograd
features. Since this is a vast topic beyond the scope of this discussion, here are some of the ways via which you can learn how to implement Automated Differentiation using AutoGrad in PyTorch:
- PyTorch AutoGrad Official Tutorial/Documentation
- Pytorch Tutorial — Gradient Descent with Autograd and Backpropagation by Python Engineer
Note : The implementation mentioned in this blog is for displaying the easiest way of differentiation — the numerical way but since it requires another function bypass call it should not be used on a large scale because a suitable time loss will occur. The below section shows the suitable implementation time differences.
Running Time of Scripts
After running the three methods separately and timing those scripts the following averaged timing data is obtained:
Over 100 Iterations :
ClassicalDifferential : 6.438255310058594e-05 seconds. NumericalDifferential : 0.00017073869705200195 seconds. TorchDifferential : 0.010568954944610597 seconds. Over 1000 Iterations :
ClassicalDifferential : 5.8611392974853516e-05 seconds. NumericalDifferential : 0.00013384795188903809 seconds. TorchDifferential : 0.010314790725708008 seconds.
Appendix: Running Time Codes —
Further Reading & Understanding
If any/all of these don’t make complete sense to you, consider reading these articles to enhance your POC:
- Brilliant — Limits of Functions
- CueMath — LHD & RHD