## Introduction

The numerical differentiation method is one among many techniques used to numerically compute derivatives. This technique is widely used due to its simplicity and accuracy. Numerical differentiation method is mostly used in computer graphics, simulation of different scientific models, numerical models for climate forecast, numerical fluid flow analysis and is especially useful where an analytical solution is not obtainable.

## Implementation

``````
(* Numerical Differentiation *)

open Owl_types

(* Functor of making numerical differentiation module of different precisions *)

module Make (A : Ndarray_Numdiff with type elt = float) = struct
type arr = A.arr

type elt = A.elt

(* global epsilon value used in numerical differentiation *)
let _eps = 0.00001

let _ep1 = 1.0 /. _eps

let _ep2 = 0.5 /. _eps

(* derivative of f : scalar -> scalar *)
let diff f x = (f (x +. _eps) -. f (x -. _eps)) *. _ep2

(* derivative of f : scalar -> scalar, return both function value and derivative *)
let diff' f x = f x, diff f x

(* second order derivative of f : scalar -> scalar *)
let diff2 f x = (f (x +. _eps) +. f (x -. _eps) -. (2. *. f x)) /. (_eps *. _eps)

(* second order derivative of f : float -> float, return both function value and derivative *)
let diff2' f x = f x, diff2 f x

(* gradient of f : vector -> scalar, return both function value and gradient *)
let n = A.numel x in
let g = A.create [| n |] (f x) in
let gg =
A.mapi
(fun i xi ->
let x' = A.copy x in
A.set x' [| i |] (A.elt_to_float xi +. _eps);
f x')
x
in
g, A.((gg - g) *\$ _ep1)

(* gradient of f : vector -> scalar *)

(* transposed jacobian of f : vector -> vector, return both function value and jacobian *)
let jacobianT' f x =
let y = f x in
let m, n = A.numel x, A.numel y in
let j = A.tile y [| m; 1 |] in
let jj = A.copy j in
for i = 0 to m - 1 do
let x' = A.copy x in
let a = A.elt_to_float (A.get x [| i |]) in
A.set x' [| i |] (a +. _eps);
let y' = A.reshape (f x') [| 1; n |] in
A.set_slice [ [ i ]; [] ] jj y'
done;
y, A.((jj - j) *\$ _ep1)

(* transposed jacobian of f : vector -> vector *)
let jacobianT f x = jacobianT' f x |> snd

(* jacobian of f : vector -> vector, return both function value and jacobian *)
let jacobian' f x =
let y, j = jacobianT' f x in
y, A.transpose ~axis:[| 1; 0 |] j

(* jacobian of f : vector -> vector *)
let jacobian f x = jacobian' f x |> snd
end

``````

The given algorithm implements the numerical differentiation method in OCaml to do the following.

• Calculate the derivative of a scalar to a scalar function
• Calculate the second-order derivative of a scalar to a scalar function
• Calculate the gradient of a vector to scalar function
• Calculate the transposed Jacobian of a vector to vector function
• Calculate the Jacobian of a vector to a vector function

These functions can also return the both function value and derivative as a tuple.

## Step-by-step Explanation

1. For the function `diff` and `diff'`, we define global values `_eps`, `_ep1`, and `_ep2` which are used in the numerical differentiation method. `_eps` is the value of Epsilon which is a small number used to approximate the derivative of a function. Optimally, it should be the square root of the machine’s precision, which is roughly 1.0e-8 in 64-bit systems. `_ep1` and `_ep2` are two values used to calculate the derivative.
2. The function `diff` accepts a scalar to scalar function and a scalar value `x` as input, and returns the derivative of that function at `x`.
• This is done by calculating the difference in function value at `x+_eps` and `x-_eps`, and then dividing it by `2*_eps`.
3. The function `diff'` is similar to `diff`, but also returns the function value along with the derivative in a tuple.
4. The function `diff2` is used to calculate the second-order derivative of a scalar to scalar function.
• Similar to `diff`, this function accepts a scalar to scalar function and a scalar value as input and then returns the second-order derivative at `x`.
5. The function `diff2'` is similar to `diff2` but also returns the function value along with the second-order derivative in a tuple.
6. The function `grad'` accepts a vector to scalar function and a vector as input, and returns the gradient of that function at the given point as a tuple.
• First, function value is evaluated at `x`. The number of elements in `x` is stored in `n`.
• We initialise a gradient vector `g` with `n` number elements and equal to the function value.
• Then we create another matrix `gg` with `n` number columns and one row for each element in `x`. This matrix represents the function values evaluated at `x+_eps` and `x-_eps`.
• For each element in `x`, the same value is copied to `x'` as in `x`. However, at the index of that element, `_eps` is added to its value.
• Evaluating the function values at each of the points `x` and `x'` gives us a matrix `gg`. By subtracting `g` from `gg` and multiplying the result with `_ep1`, we eventually obtain the gradient.
7. The function `grad` is similar to `grad'`, but it only returns the calculated gradient.
8. The function `jacobianT'` is used to calculate the transposed Jacobian matrix of a vector to vector function.
• First, the value of the function `f` is evaluated at `x`.
• `j` is created as a `m x n` matrix where m is the number of elements in the vector `x` and n is the number of elements in the function evaluated at `x`.
• We then copy `j` to `jj` so that we have a copy of the matrix j to which we will apply the finite difference method.
• For each element in `x`, we create a copy of `x` and add `_eps` to the element at ith position to get `x'`.
• The function is evaluated at `x'` and reshaped into an array of size `1 x n`. Then it is stored at the ith row of `jj`.
• We subtract the matrix `j` from the resulting matrix `jj` and multiply the result with `_ep1` to get the transposed Jacobian of f.
9. The function `jacobianT` is similar to `jacobianT'`, but only returns the transpose Jacobian of the function.
10. The function `jacobian'` is used to calculate the Jacobian matrix of a vector to vector function.
• First, we calculate the value of the function `f` at `x` and store it in variable `y`.
• We then calculate the transpose jacobian of `f` at `x` and store it in variable `j`.
• We return a tuple of `y` and the transposition of `j`.
11. The function `jacobian` is similar to `jacobian'`, but only returns the Jacobian of the function.

## Complexity Analysis

The time complexity of the given algorithm is dependent on the input function and the dimension of the input function parameters. The gradient and Jacobian functions have time complexity of O(n^2) where n is the dimension of the input, since they involve multiple evaluations of the input function. The functions `diff`, `diff'`, `diff2` and `diff2'` has time complexity of O(1), since they only require two evaluations of the input function. Due to the epsilon value used to approximate the derivative, these functions have a reduced accuracy compared to mathematical differentiation and may require a smaller value of epsilon to compute with better accuracy.