static Tutorial 02: Computing derivatives.

Contributors:
    Mauricio Aristizabal
        University of Texas at San Antonio
        Universidad EAFIT, Colombia.

Date of creation:  Jun 20 2024
Last modification: Jun 20 2024

Introduction

The main purpose of this document is to show how to compute derivatives of functions in python using the static ‘Order truncated Imaginary’ (OTI) module in pyoti. Multiple examples are given in order to show and understand the capabilities of the library, as well as its basic advantages. For basic usage of the static module, see the static Tutorial 01.

For this, first import the libraries required.

import pyoti.static.onumm2n2 as oti  # Import the library for m=2 and n=3
import numpy as np

Derivatives of single-variable functions.

In general, OTI numbers have a strong application in finding derivatives of functions. In this sense let’s see some applications to obtain derivatives.

Consider a \(C^n\) function \(f: \mathbb{R} \rightarrow \mathbb{R}\), and it’s Taylor series:

\[f(z) = \sum_{p=0}^{\infty} \frac{1}{p!}\frac{\partial^p f}{\partial x^p}(a)\;(z-a)^p\]

Let´s take, for convenience, \(z=x_0+\epsilon_1\) and \(a=x_0\). Then the following is obtained:

\[f(x_0+\epsilon_1 ) = f(x_0) + \displaystyle \frac{\partial f}{\partial x }(x_0)\;{\epsilon_1 } + \displaystyle\frac{1}{2!}\frac{\partial^2 f}{\partial x^2}(x_0)\;{\epsilon_1^2} + \displaystyle\frac{1}{3!}\frac{\partial^3 f}{\partial x^3}(x_0)\;{\epsilon_1^3} + ...\]

Therefore, the result of evaluating a function at an OTI number in the form \(x_0+\epsilon_1\), is an OTI number that contains the derivatives of the function along its imaginary directions. Also note that the previous expansion will go until the truncation order set to the number \(x_0+\epsilon_1\). Therefore, the truncation order \(n\) is set to the maximum order of derivative required to be computed. Also, typically, the notation for a perturbed variable is \(x^* = x_0+\epsilon_1\)

Notice two things: (i) that the coefficient along the real part contains the function evaluated at \(x_0\); and (ii) that one gets the \(p\)’th order derivative evaluated at \(x_0\) multiplied by the factor \(p!\) in the coefficient along the \(p\)’th order \(\epsilon_1\) direction.

Therefore, to retrieve a derivative, one must extract the imaginary coefficient along the \(\epsilon_1^p\) direction and multiply (remove) the multiplication factor, i.e.

\[\begin{split}\begin{array}{c} \displaystyle\frac{\partial f}{\partial x}(x_0) = \mbox{Im}_{\epsilon_1}\left[f(x_0+\epsilon_1)\right] \\ \displaystyle\frac{\partial^2 f}{\partial x^2}(x_0) = 2! \mbox{Im}_{\epsilon_1^2}\left[f(x_0+\epsilon_1)\right] \\ \displaystyle\frac{\partial^3 f}{\partial x^3}(x_0) = 3! \mbox{Im}_{\epsilon_1^3}\left[f(x_0+\epsilon_1)\right] \\ \end{array}\end{split}\]

and in general

\[\displaystyle\frac{\partial^p f}{\partial x^p}(x_0) = p! \mbox{Im}_{\epsilon_1^p}\left[f(x_0+\epsilon_1)\right]\]

To apply this in a real example let’s consider the following function:

Example: Polynomial function

Consider the following function:

\[f(x)= 10x^5-3x^3+2x-10\]

We will compute derivatives up to 5th order.

def funct_1(x):
    return 10.0*x**5 -3.0*x**3+2.0*x-10.0

The analytical defivatives of the funciton are:

\[\begin{split}\begin{array}{rcl} \displaystyle\frac{\partial f}{\partial x }(x) &=& 50x^4-9x^2+2, \\ \displaystyle\frac{\partial^2 f}{\partial x^2}(x) &=& 200x^3-18x, \\ \displaystyle\frac{\partial^3 f}{\partial x^3}(x) &=& 600x^2-18, \\ \end{array}\end{split}\]

Let’s use this simple exampel and use OTIs to compute up-to 5th order derivatives of this polynomial.

First, apply the perturbation to the variable \(x\) along \(\epsilon_1\), \(x^* = 2.5+\epsilon_1\). Since want to obtain higher order terms, we need to redefine the order of the number. Therefore we will redefine it to order 5 (because we want to get up to the fifth derivative).

x = 2.5
x_star= x + oti.e(1)

print(f'Variable x:       {x_star}')
print(f'Truncation order: {x_star.order}')

Variable x:       2.5000 + 1.0000 * e([1]) + 0.0000 * e([2]) + 0.0000 * e([[1,2]]) + 0.0000 * e([1,2]) + 0.0000 * e([[2,2]])
Truncation order: 2

Now let’s evaluate the function at \(x^*\)

result_1 = funct_1(x_star)
print(result_1)

924.6875 + 1898.8750 * e([1]) + 0.0000 * e([2]) + 1540.0000 * e([[1,2]]) + 0.0000 * e([1,2]) + 0.0000 * e([[2,2]])

Observe that the first two values correspond explicitly to the function and first derivative evaluated at \(x_0=2.5\). Only the values along \(\epsilon_1\) or its multipls are non-zero.

In order to retreive the values of the derivatives, the solution must take into account the following:

\[\mbox{result} = a_r + a_{\epsilon_1}\epsilon_1 + a_{\epsilon_2}\epsilon_2 + a_{\epsilon_1^2}\epsilon_1^2 + + a_{\epsilon_1\epsilon_2}\epsilon_1\epsilon_2 + a_{\epsilon_2^2}\epsilon_2^2\]

and considering that

\[f(x_0+\epsilon_1,y_0+\epsilon_2) = f(x_0,y_0) + \displaystyle \frac{\partial f}{\partial x }(x_0,y_0){\epsilon_1 } + \displaystyle \frac{\partial f}{\partial y }(x_0,y_0){\epsilon_2 } + \displaystyle\frac{1}{2!}\frac{\partial^2 f}{\partial x^2}(x_0,y_0){\epsilon_1^2} + \displaystyle \frac{\partial^2 f}{\partial x\partial y}(x_0,y_0){\epsilon_1\epsilon_2} + \displaystyle\frac{1}{2!}\frac{\partial^2 f}{\partial y^2}(x_0,y_0){\epsilon_2^2} + ...\]

then to retreive the derivatives, we most do the following:

\[\begin{split}\begin{array}{rcl} f(x_0) &=& a_r \\ \displaystyle \frac{\partial f}{\partial x }(x_0,y_0)) &=& a_{\epsilon_1} \\ \displaystyle \frac{\partial f}{\partial y }(x_0,y_0) &=& a_{\epsilon_2} \\ \displaystyle\frac{\partial^2 f}{\partial x^2}(x_0,y_0)&=& a_{\epsilon_1^2}\;2! \\ \displaystyle\frac{\partial^2 f}{\partial x\partial y}(x_0,y_0)&=& a_{\epsilon_1\epsilon_2} \\ \displaystyle\frac{\partial^2 f}{\partial y^2}(x_0,y_0)&=& a_{\epsilon_2^2}\;2! \\ \end{array}\end{split}\]

As a result, in order to evaluate the result, consider the following:

f       = result_1.real
f_x     = result_1.get_im([1])
f_xx    = result_1.get_im([[1,2]]) * 2*1

An easier way to do this is to use the method get_deriv(dirArray). This automatically applies the factorial factor to obtain the derivative.

print("f_xx:")
print(result_1.get_im([1,1]) * 2*1, # Using get_im and the factorial
      result_1.get_deriv([1,1])     # Using just get_deriv.
     )

f_xx:
3080.0 3080.0

# Getting the other derivatives
f_xxx  = result_1.get_deriv([[1,3]])

# Print results:
print(f"Derivative |   pyoti  | analytical|")
print(f"-----------|----------|-----------|")
print(f"f          | {f:8g} | {(10*x**5-3*x**3+2*x-10):9g} |")
print(f"f_x        | {f_x:8g} | {(50*x**4-9*x**2+2):9g} |")
print(f"f_xx       | {f_xx:8g} | {(200*x**3-18*x):9g} |")    # ")#
print(f"f_xxx      | {f_xxx:8g} | {(600*x**2-18):9g} |")

Derivative |   pyoti  | analytical|
-----------|----------|-----------|
f          |  924.688 |   924.688 |
f_x        |  1898.88 |   1898.88 |
f_xx       |     3080 |      3080 |
f_xxx      |        0 |      3732 |

Derivatives of multi variable functions

A similar approach as in the single variable case applies to multivariable functions. If one considers a function \(f: \mathbb{R}^m \rightarrow \mathbb{R}\), then all its derivatives can be obtained by perturbing the input parameters along independent imaginary directions, as follows:

\[\begin{split}\begin{array}{rcl} x_1^* &=& x_1 + \epsilon_1\\ x_2^* &=& x_2 + \epsilon_2\\ &\vdots& \\ x_m^* &=& x_m + \epsilon_m\\ \end{array}\end{split}\]

Then, the function is evaluated using the OTI function \(f(x_1^*,x_2^*,...,x_m^*)\). The derivatives are obtained as follows

\[\frac{ \partial^{p}f }{ \partial x_1^{p_1} \partial x_2^{p_2} ... \partial x_m^{p_m} } = p_1!p_2! \cdots p_m! \; \mbox{Im}_{\epsilon_1^{p_1}\epsilon_2^{p_2}\cdots\epsilon_m^{p_m}}\left[f(x_1^*,x_2^*,...,x_m^*)\right]\]

where \(p=\sum_{i=1}^m p_i\)

For example, a two variable function up to order 2 derivatives gives:

\[\begin{split}f(x_{1_0}+\epsilon_1,x_{2_0}+\epsilon_2) = f(x_{1_0},x_{2_0}) + \frac{\partial f }{ \partial x_1}(x_{1_0},x_{2_0}) \ \epsilon_1 + \frac{\partial f }{ \partial x_2}(x_{1_0},x_{2_0}) \ \epsilon_2 + \\ \frac{1}{2!}\frac{\partial^2 f }{ \partial x_1^2}(x_{1_0},x_{2_0}) \ \epsilon_1^2 + \frac{\partial^2 f }{ \partial x_1\partial x_2}(x_{1_0},x_{2_0}) \ \epsilon_2\epsilon_1 + \frac{1}{2!}\frac{\partial^2 f }{ \partial x_2^2}(x_{1_0},x_{2_0}) \ \epsilon_2^2 + \cdots\end{split}\]

Example

To apply this to a simple example, let’s consider the function

\[f(x,y) = cos( (x+2y-1)^2 )\]

Let’s define the function in the python environment:

def funct_2(x,y):
    return oti.cos( ( x + 2*y - 1 ) ** 2 )

Analytically, let’s define the derivatives of the function:

\[\begin{split}f_x(x,y) = -2(x+2y-1)\sin((x+2y-1)^2) , \\ f_y(x,y) = -4(x+2y-1)sin((x+2y-1)^2)\end{split}\]

Now, how do we retreive the function’s first derivatives by just evaluating the function with OTI numbers? It’s easy!

According to the taylor series expansion of the function, we can retreive the function derivatives by evaluating the function at the point we want, say \((x_0,y_0)\) with a perturbation along OTI directions. In this case, the perturbation will occur as follows: \((x_0+\epsilon_1,y_0+\epsilon_2)\)

To get all the derivatives in a single evaluation of the function derivatives one most only evaluate the function into the OTI numbers. In this way, we define the point \((x_0=2,y_0=3)\) and we get:

x=2.0 + oti.e(1)
y=3.0 + oti.e(2)

result_2 = funct_2(x,y)
print(result_2)

0.3006 + 13.3525 * e([1]) + 26.7051 * e([2]) - 28.5043 * e([[1,2]]) - 114.0173 * e([1,2]) - 114.0173 * e([[2,2]])

All derivatives are now computed and lay in the coefficients of teh OTI number.

Comparing against the analytical derivatives, the following results are obtained:

f_an = np.cos( ( x.real + 2*y.real - 1 ) ** 2 )
print(f' f(x,y): analytical: {f_an}')
print(f'              pyoti: {result_2.real}')
print(f'              Error: {f_an-result_2.real}')

 f(x,y): analytical: 0.30059254374363714
              pyoti: 0.3005925437436371
              Error: 5.551115123125783e-17

The results of the first-order derivatives are as follows:

f_x_an = -2*( x.real + 2*y.real - 1 ) * np.sin( ( x.real + 2*y.real - 1 ) ** 2 )
print(f' f_x(x,y): analytical: {f_x_an}')
print(f'                pyoti: {result_2.get_deriv([1])}')
print(f'                Error: {f_x_an-result_2.get_deriv([1])}')

 f_x(x,y): analytical: 13.352537138632607
                pyoti: 13.352537138632606
                Error: 1.7763568394002505e-15

f_y_an = -4*( x.real + 2*y.real - 1 ) * np.sin( ( x.real + 2*y.real - 1 ) ** 2 )
print(f' f_y(x,y): analytical: {f_y_an}')
print(f'                pyoti: {result_2.get_deriv([2])}')
print(f'                Error: {f_y_an-result_2.get_deriv([2])}')

 f_y(x,y): analytical: 26.705074277265215
                pyoti: 26.70507427726521
                Error: 3.552713678800501e-15

The seccond derivatives of the function are:

\[\begin{split}f_{xx}(x,y) = -2\sin((x+2y-1)^2)-4(x+2y-1)^2\cos((x+2y-1)^2), \\ f_{xy}(x,y) = -4\sin((x+2y-1)^2)-8(x+2y-1)^2\cos((x+2y-1)^2), \\ f_{yy}(x,y) = -8\sin((x+2y-1)^2)-16(x+2y-1)^2\sin((x+2y-1)^2)\end{split}\]

f_xx_an = -2*oti.sin((x.real+2*y.real-1)**2)-4*(x.real+2*y.real-1)**2*np.cos((x.real+2*y.real-1)**2)
print(f' f_xx(x,y): analytical: {f_xx_an}')
print(f'                 pyoti: {result_2.get_deriv([1,1])}')
print(f'                 Error: {f_xx_an-result_2.get_deriv([1,1])}')

 f_xx(x,y): analytical: -57.008633268233936
                 pyoti: -57.00863326823392
                 Error: -1.4210854715202004e-14

f_xy_an = -4*oti.sin((x.real+2*y.real-1)**2)-8*(x.real+2*y.real-1)**2*np.cos((x.real+2*y.real-1)**2)
print(f' f_xy(x,y): analytical: {f_xy_an}')
print(f'                 pyoti: {result_2.get_deriv([1,2])}')
print(f'                 Error: {f_xy_an-result_2.get_deriv([1,2])}')

 f_xy(x,y): analytical: -114.01726653646787
                 pyoti: -114.01726653646784
                 Error: -2.842170943040401e-14

f_yy_an = -8*oti.sin((x.real+2*y.real-1)**2)-16*(x.real+2*y.real-1)**2*np.cos((x.real+2*y.real-1)**2)
print(f' f_yy(x,y): analytical: {f_yy_an}')
print(f'                 pyoti: {result_2.get_deriv([2,2])}')
print(f'                 Error: {f_yy_an-result_2.get_deriv([2,2])}')

 f_yy(x,y): analytical: -228.03453307293574
                 pyoti: -228.0345330729357
                 Error: -5.684341886080802e-14