LChrisman/Computing Derivatives for MetaLog

This page outlines an algorithm for computing the [math]\displaystyle{ i^{th} }[/math] derivative of a Keelin MetaLog.

From Eq(6) in the Baucells, Chrisman, Keelin paper draft:

[math]\displaystyle{ M^{(i)}(y) = \mu^{(i)}(y) + s^{(i)} \ell(y) + \sum_{j=0}^{i-1} s^{(j)}(y) \ell^{(i-j)}(y) }[/math]

thus in

[math]\displaystyle{ y^i (1-y)^i M^{(i)}(y) = y^i (1-y)^i \mu^{(i)}(y) + y^i (1-y)^i s^{(i)} \ell(y) + y^i (1-y)^i \sum_{j=0}^{i-1} s^{(j)}(y) \ell^{(i-j)}(y) }[/math]

the first and third term are polynomials, and the middle term is a polynomial times a logit. Hence, this can always be written in the form:

[math]\displaystyle{ y^i (1-y)^i M^{(i)}(y) = poly_1(y) + poly_2(y) \ell(y) }[/math]

where

[math]\displaystyle{ poly_2(y) = y^i (1-y)^i s^{(i)} }[/math]

[math]\displaystyle{ poly_1(y) = y^i (1-y)^i \mu^{(i)}(y) + T(y) }[/math]

[math]\displaystyle{ T(y) = y^i (1-y)^i \sum_{j=0}^{i-1} s^{(j)}(y) \ell^{(i-j)}(y) }[/math]

When [math]\displaystyle{ i\ge\lfloor (k+1)/2 \rfloor }[/math], then [math]\displaystyle{ poly_2(y)=0 }[/math] so that [math]\displaystyle{ y^i (1-y)^i M^{(i)}(y) }[/math] is a simple polynomial.

Change of variable

[math]\displaystyle{ u = y-0.5 }[/math]

[math]\displaystyle{ T(u) = T(y-0.5) }[/math]

i^th derivative finding algorithm

The algorithm uses arithmetic operations on a polynomial-coefficients data type. Here [math]\displaystyle{ T(u) }[/math], [math]\displaystyle{ s^{(j)} }[/math], [math]\displaystyle{ R(u) }[/math] and [math]\displaystyle{ N(u) }[/math] are all polynomials. The order of the derivative, [math]\displaystyle{ i }[/math] is given.

1. For each [math]\displaystyle{ j=0 \mbox{ to } i }[/math]
   • Compute [math]\displaystyle{ C_j = \left( \begin{array}{c}i \\ j\end{array}\right) ( i-j-1)! }[/math]
   • Compute [math]\displaystyle{ s^{(j)}(u) }[/math]
   • Compute [math]\displaystyle{ R(u) = (0.5+u)^j (0.5-u)^j }[/math]
   • Compute [math]\displaystyle{ N(u) = (0.5+u)^{i-j} - (-1)^{i-j} ( 0.5-u)^{i-j} }[/math]
   • Compute [math]\displaystyle{ T_j(u) = C_j \cdot s^{(j)}(u) \cdot R(u) \cdot N(u) }[/math]
2. [math]\displaystyle{ T(u) = \sum_{j=0}^i T_j(u) }[/math]

Zero finding

To find the zeros of [math]\displaystyle{ M^{(i)}(y) }[/math], first find the zeros of [math]\displaystyle{ M^{(i+1)}(y) }[/math], which are the minima and maxima of [math]\displaystyle{ M^{(i)}(y) }[/math]. We should add to this list of extrema y=0 and y=1. Each adjacent extrema with opposite sign will bound one zero of [math]\displaystyle{ M^{(i)}(y) }[/math]. Thus, we have a finite set of intervals, each containing a single zero.

The base case for the recursion is [math]\displaystyle{ i=\lfloor (k+1)/2 \rfloor }[/math], in which case [math]\displaystyle{ y^i (1-y)^i M'(y) = y^i (1-y)^i \mu^{(i)} + T(y) }[/math] is a polynomial, for which a standard polynomial root finder is used. This has the same (interior) roots as [math]\displaystyle{ M'(y) }[/math].

To find the zero in each interval, simply conduct a binary search within each interval using [math]\displaystyle{ y^i (1-y)^i M^{(i)}(y) = y^i (1-y)^i s^{(i)} + T(y) + y^i (1-y)^i s^{(i)} \ell(y) }[/math]

It is possible to speed the search up within each interval with a bounded Newton algorithm. For this, we need the derivative

[math]\displaystyle{ {{\partial y^i (1-y)^i M^{(i)}(y)} \over {\partial y}} = {{\partial y^i (1-y)^i}\over {\partial y}} M^{(i)}(y) + y^i (1-y)^i M^{(i+1)}(y) }[/math]

where [math]\displaystyle{ y^i (1-y)^i M^{(i+1)}(y) }[/math] is obtained from the previous recursion step (at [math]\displaystyle{ i+1 }[/math]).