## 2 Background
### Arbitrary-Order Differentiation Using HYPAD  
HYPAD arises from the Complex Taylor Series Expansion (CTSE), a first-order numerical differentiation technique similar in spirit and concept to finite differencing but with significant numerical advantages. CTSE uses the orthogonality of the real and imaginary axes of the complex plane to calculate derivatives. Unlike finite differencing, CTSE does not require the difference between two analyses; hence, the specification of step size is a non-issue. Instead, the derivatives are computed using a perturbation along the imaginary axis. Dual numbers that are a subset of the generalized complex numbers, also known as hypercomplex [67], can be used for the non-real perturbation in an equivalent fashion to complex [68]. With the use of dual numbers, both truncation and subtractive cancellation errors are eliminated in contrast to traditional complex where only the truncation error is eliminated [49]. An important feature of such is that any perturbation magnitude can be used without losing accuracy. For instance, to obtain the first-order derivatives for a function \(f(\theta)\), the variable \(\theta\) is perturbed to become \(\theta^{*}=\theta+h\epsilon_{1}\), where \(\epsilon_{1}\) is the imaginary component of the dual number, \(h\) denotes the step size, and \(\theta^{*}\) is the dual representation of order one of \(\theta\). As \(f\) becomes a dual function, it can be expanded using CTSE as follows:  
\[f(\theta+h\epsilon_{1})=f(\theta)+\frac{h\epsilon_{1}}{1!}\frac{df(\theta)}{ d\theta}+\frac{(h\epsilon_{1})^{2}}{2!}\frac{d^{2}f(\theta)}{d\theta^{2}}+ \frac{(h\epsilon_{1})^{3}}{3!}\frac{d^{3}f(\theta)}{d\theta^{3}}+\ldots \tag{5}\]  
Using the properties of dual numbers, i.e., \(\epsilon_{1}^{q}=0\) for \(q\geq 2\), and \(\epsilon\
eq 0\), the terms \(\frac{(h\epsilon_{1})^{2}}{2! }\frac{d^{2}f(\theta)}{d\theta^{2}},\frac{(h \epsilon_{1})^{3}}{3! }\frac{d^{3}f(\theta)}{d\theta^{3}}\) and the rest of the higher order terms are eliminated. Therefore, the first-order derivative can be exactly estimated by taking the dual imaginary part of both sides from Equation (5) as:  
\[\frac{df(\theta)}{d\theta}=\frac{Im_{1}[f(\theta+h\epsilon_{1})]}{h} \tag{6}\]  
where \(Im_{1}[]\) denotes the dual imaginary component of the dual number corresponding to the axis \(\epsilon_{1}\). Using this approach, the error in the derivative approximation can be reduced to machine precision because no subtraction cancellation error is present. In addition, the method has no truncation error; therefore, there is no dependence on the parameter \(h\)[69]. In general, a step size of \(h=1\) is used for simplicity. Therefore, \(h\) is omitted in the expressions presented hereafter. Extending CTSE to hypercomplex mathematics gives rise to the hypercomplex automatic differentiation (HYPAD), which mathematically allows for the computing of machine precision higher-order and mixed derivatives up to any order [48; 49; 70]. Conceptually the CTSE and HYPAD are based on the same mathematical fundamentals. The only difference resides in the use of multidual numbers in HYPAD to compute the hypercomplex Taylor series expansion of a function and obtain arbitrary-order sensitivities. In a nutshell, a multidual number consists of multiple independent dual imaginary axes, e.g., \(\epsilon_{1},\epsilon_{2},\ldots,\epsilon_{m}\). These multiple imaginary axes are independent and orthogonal, and high-order and mixed partial derivatives are obtained by perturbing the parameters of interest along the different independent imaginary axes. As such, the order of the derivatives, \(m\), defines the order of the multidual number required to evaluate the function. The total number of imaginary axes for a \(m^{th}\) order multidual number is given by \(p=2^{m}-1\). For instance, to obtain second-order derivatives, \(m=2\) and \(p=3\); as such, three imaginary axes (e.g., \(\epsilon_{1},\epsilon_{2}\), and \(\epsilon_{12}\)) are needed, where two are independent and one is a mixed axis.