Automatic differentiation using dual numbers

[Albert Chan]

(06-22-2022 10:36 PM)Thomas Klemm Wrote:  Halley's Method

It's straightforward to implement Halley's method:

\(
\begin{align}
x_{{n+1}}=x_{n}-{\frac{2f(x_{n})f'(x_{n})}{2{[f'(x_{n})]}^{2}-f(x_{n})f''(x_{n})}}
\end{align}
\)

However we rather use this formula:

\(
\begin{align}
x_{n+1}
&=x_{n}-{\frac {f(x_{n})}{f'(x_{n})-{\frac {f(x_{n})}{f'(x_{n})}}{\frac {f''(x_{n})}{2}}}} \\
&=x_{n}-{\frac {f(x_{n})}{f'(x_{n})}}\left[1-{\frac {f(x_{n})}{f'(x_{n})}}\cdot {\frac {f''(x_{n})}{2f'(x_{n})}}\right]^{-1} \\
\end{align}
\)

Is correction factor to Newton's always positive ?
How do we ensure Halley's method is safe to use ?

I've seen some page (search for householder) flip the correction: 1/(1-ε) ≈ 1+ε
Correction is less aggressive, but safer (if f'≠0, no division by zero problem)

\( \displaystyle
x_{n+1}
=x_{n}-{\frac {f(x_{n})}{f'(x_{n})}}\left[1+{\frac {f(x_{n})}{f'(x_{n})}}\cdot {\frac {f''(x_{n})}{2f'(x_{n})}}\right]
\)

[Albert Chan]

An example that c = (1 - 1/2 * (f/f') * (f''/f')) can turn negative.

f = x^4 - 30*x^3 +63*x^2 - 12*x + 33      // real roots ≈ [2.31090, 27.7432]
f' = 4*x^3 - 90*x^2 + 126*x - 12              // real roots ≈ [0.102744, 1.38994, 21.0073]
f'' = 12*x^2 -180*x + 126                         // real roots ≈ [0.736125, 14.2639]

m = f'^2 - 1/2*f*f''                                    // real roots ≈ [-0.34128, 0.484147]
c = m / f'^2

c = ±∞ when guess = f' roots (Newton's correction also ±∞)
c = 0 when guess = m roots (Halley's method hit by divide-by-zero)
c < 0 when guess within m roots (iterate opposite direction of Newton's)

Comment:

Above example is easier to solve for roots of f(1/x), then flip the results.

g = 33*x^4 - 12*x^3 + 63*x^2 - 30*x + 1
g' = 132*x^3 - 36*x^2 + 126*x - 30
g'' = 396*x^2 - 72*x + 126

g'', with negative quadratic discriminant, have no root
g'' > 0 ⇒ g is concave up.

[Albert Chan]

If guess is close to root, we can consider f as the independent variable.

x = x0 + (dx/df)(f=f0) * (f-f0)/1! + (d²x/df²)(f=f0) * (f-f0)²/2! + ...

dx/df = 1 / (df/dx) = 1/f'
d²x/df² = d/df (dx/df) = d/dx (1/f') * dx/df = (-f''/f'^2) / f' = -f''/f'^3

Solving root is simplify x at f=0, with iteration formula:

x ← x + (1/f')*(-f) + (-f''/f'^3)*(f^2/2) + ... = x - f/f' - (f^2*f'')/(2*f'^3) + ...

Drop ... = O(f^3), we have:

(06-25-2022 11:05 AM)Albert Chan Wrote:  I've seen some page (search for householder) flip the correction: 1/(1-ε) ≈ 1+ε
Correction is less aggressive, but safer (if f'≠0, no division by zero problem)

\( \displaystyle
x_{n+1}
=x_{n}-{\frac {f(x_{n})}{f'(x_{n})}}\left[1+{\frac {f(x_{n})}{f'(x_{n})}}\cdot {\frac {f''(x_{n})}{2f'(x_{n})}}\right]
\)

[Ángel Martin]

(06-26-2022 10:55 AM)Albert Chan Wrote:  

"Concave and Convex" works better for me.

[Ángel Martin]

(06-21-2022 02:19 AM)Thomas Klemm Wrote:  

(06-20-2022 09:50 PM)Albert Chan Wrote:  … forget about ε … ε is a bad idea

Au contraire!
It's just that you misinterpreted the meaning of the coefficient of \(\varepsilon^2\).
Based on the Taylor series we have:

\(
\begin{align}
f(x + \varepsilon) = f(x) + {f}'(x) \cdot \varepsilon + \frac{1}{2} {f}''(x) \cdot \varepsilon^2
\end{align}
\)

Given your example:

Quote:{f, f', f''} = {1,2,3}
{g,g',g''} = {3,1,4}

… we have to calculate:

\(
\begin{align}
\left(1 + 2 \varepsilon + \frac{3}{2} \varepsilon^2 \right)\left(3 + \varepsilon + 2 \varepsilon^2 \right)
&= 3 + 7 \varepsilon + \frac{17}{2} \varepsilon^2 + \frac{11}{2} \varepsilon^3 + 3 \varepsilon^4 \\
&= 3 + 7 \varepsilon + \frac{17}{2} \varepsilon^2 \\
\end{align}
\)

… since \(\varepsilon^3 = 0\).
With this interpretation we get the correct result.

This can be extended in an analogous way for higher-order derivatives.

Sorry for coming this late to this, but I'm not sure I follow your point about higher order derivatives (I'm sure it's my own fault). By definition \(\varepsilon^2 = 0\)., right?

[Thomas Klemm]

(07-07-2022 07:40 AM)Ángel Martin Wrote:  By definition \(\varepsilon^2 = 0\)., right?

If we want to extend the idea of dual numbers for higher order derivatives we can define \(\varepsilon^n = 0\).
Slightly sloppy, we ignore the terms of \(\varepsilon^n\) and higher.

[ijabbott]

(07-07-2022 07:35 AM)Ángel Martin Wrote:  

(06-26-2022 10:55 AM)Albert Chan Wrote:  [ "Concave up / Concave down" illustration was here ]

"Concave and Convex" works better for me.

Then you need to remember the viewpoint.