sparse Tutorial 02: Computing arbitrary-order derivatives.

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

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

Introduction

The main purpose of this document is to show how to compute derivatives of functions in python using the sparse ‘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 this, first import the libraries required.

import pyoti.sparse as oti  # Import the library.
import numpy as np

# Set to print all components
oti.set_printoptions(terms_print=-1)

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 what we get is the following:

\[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, \\ \displaystyle\frac{\partial^4 f}{\partial x^4}(x) &=& 1200x, \\ \displaystyle\frac{\partial^5 f}{\partial x^5}(x) &=& 1200, \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],order=5)

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

Variable x:       2.5 + 1 * e([1])
Truncation order: 5

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

result_1 = funct_1(x_star)
print(result_1)

924.688 + 1898.88 * e([1]) + 1540 * e([[1,2]]) + 622 * e([[1,3]]) + 125 * e([[1,4]]) + 10 * e([[1,5]])

Observe that the first two values correspond explicitly to the function and first derivative evaluated at \(x_0=2.5\). In order to retreive the values of the other derivatives, the solution must take into account the following:

if our number results in something like

\[\mbox{result} = a_0 + a_1\epsilon_1 + a_2\epsilon_1^2 + a_3\epsilon_1^3 + a_4\epsilon_1^4 + a_5\epsilon_1^5\]

and considering that

\[f(x_0+\epsilon) = 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} + ...\]

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

\[\begin{split}\begin{array}{rcl} f(x_0) &=& a_0 \\ \displaystyle \frac{\partial f}{\partial x }(x_0) &=& a_1 \\ \displaystyle\frac{\partial^2 f}{\partial x^2}(x_0)&=&a_2*2! \\ \displaystyle\frac{\partial^3 f}{\partial x^3}(x_0)&=&a_3*3! \\ \displaystyle\frac{\partial^4 f}{\partial x^4}(x_0)&=&a_4*4! \\ \displaystyle\frac{\partial^5 f}{\partial x^5}(x_0)&=&a_5*5! \\ \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
f_xxx   = result_1.get_im([[1,3]]) * 3*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_xxx:")
print(result_1.get_im([[1,3]]) * 3*2*1, # Using get_im and the factorial
      result_1.get_deriv([[1,3]])       # Using just get_deriv.
     )

f_xxx:
3732.0 3732.0

# Getting the other derivatives
f_xxxx  = result_1.get_deriv([[1,4]])
f_xxxxx = result_1.get_deriv([[1,5]])

# 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} |")
print(f"f_xxxx     | {f_xxxx:8g} | {(1200*x):9g} |")
print(f"f_xxxxx    | {f_xxxxx:8g} | {1200:9g} |")

Derivative |   pyoti  | analytical|
-----------|----------|-----------|
f          |  924.688 |   924.688 |
f_x        |  1898.88 |   1898.88 |
f_xx       |     3080 |      3080 |
f_xxx      |     3732 |      3732 |
f_xxxx     |     3000 |      3000 |
f_xxxxx    |     1200 |      1200 |

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) = 4x^5y^4\]

Let’s define the function in the python environment:

def funct_2(x,y):
    return 4*x**5*y**4

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

\[f_x(x,y) = 20x^4y^4, \ \ f_y(x,y) = 16x^5y^3\]

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 first derivatives we most only evaluate the function into OTI numbers of order one. In this way, we define the point \((x_0=2,y_0=3)\) and we get:

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

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

10368 + 25920 * e([1]) + 13824 * e([2]) + 25920 * e([[1,2]]) + 34560 * e([1,2]) + 6912 * e([[2,2]])

Analitically, the results of evaluating all function and its derivatives are:

\[f(2,3) = 10368\]

\[f_x(2,3) = 25920\]

\[f_y(2,3) = 13824\]

Now, interestingly, the coefficients that we obtain are the following: 2048.0 on the real part, wich represents the function evaluated at \((2,3)\), the coefficient of \(\epsilon_1\) is the first derivative with respect to \(x\) (as is) and the coefficient of \(\epsilon_2\) is the first derivative with respect to \(y\) function evaluated.

The seccond derivatives of the function are:

\[\begin{split}f_{xx}(x,y) = 80x^3y^4, \\ f_{xy}(x,y) = 80x^4y^3, \\ f_{yy}(x,y) = 48x^5y^2\end{split}\]

Wich evaluated at \((2,3)\) gives

\[\begin{split}f_{xx}(2,3) = 51840, \\ f_{xy}(2,3) = 34560, \\ f_{yy}(2,3) = 13824\end{split}\]

Lets now change the order of both numbers, and see what we get:

# One can increase dynamically the order of an OTI number by
# adding a zero-coefficient OTI numbers with the new desired
# truncation order.
x = x + 0*oti.e(1,order=3)
y = y + 0*oti.e(1,order=3)

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

10368 + 25920 * e([1]) + 13824 * e([2]) + 25920 * e([[1,2]]) + 34560 * e([1,2]) + 6912 * e([[2,2]]) + 12960 * e([[1,3]]) + 34560 * e([[1,2],2]) + 17280 * e([1,[2,2]]) + 1536 * e([[2,3]])

We can extract, based on the Taylor series, the different terms. However, in order to get an easy form to extract them, the method get_deriv(dirArray) helps in retreiving the derivatives according to the specifyed direction. It works as follows:

print(f"f(x_0, y_0)    :{result_2_order3.get_deriv([0]):8g}")
print()
print(f"f_x(x_0, y_0)  :{result_2_order3.get_deriv([1]):8g}")
print(f"f_y(x_0, y_0)  :{result_2_order3.get_deriv([2]):8g}")
print()
print(f"f_xx(x_0, y_0) :{result_2_order3.get_deriv([[1,2]]):8g}")
print(f"f_xy(x_0, y_0) :{result_2_order3.get_deriv([1,2]):8g}")
print(f"f_yy(x_0, y_0) :{result_2_order3.get_deriv([[2,2]]):8g}")
print()
print(f"f_xxx(x_0, y_0):{result_2_order3.get_deriv([1,1,1]):8g}")
print(f"f_xxy(x_0, y_0):{result_2_order3.get_deriv([1,1,2]):8g}")
print(f"f_xyy(x_0, y_0):{result_2_order3.get_deriv([1,2,2]):8g}")
print(f"f_yyy(x_0, y_0):{result_2_order3.get_deriv([2,2,2]):8g}")

f(x_0, y_0)    :   10368

f_x(x_0, y_0)  :   25920
f_y(x_0, y_0)  :   13824

f_xx(x_0, y_0) :   51840
f_xy(x_0, y_0) :   34560
f_yy(x_0, y_0) :   13824

f_xxx(x_0, y_0):   77760
f_xxy(x_0, y_0):   69120
f_xyy(x_0, y_0):   34560
f_yyy(x_0, y_0):    9216