Automatic differentiation using dual numbers

[Albert Chan]

(05-17-2022 01:31 PM)Albert Chan Wrote:  lua> n, pv, pmt, fv = 10, 50, -30, 100
lua> g = loan_rate2(n,pv,pmt,fv)
lua> g()
-0.2844446981069045        0.015555301893095532
lua> g()
-0.28443599888025756      8.699226646944535e-006
lua> g()
-0.28443599888025595      1.6279605736411871e-015

Perhaps we can extend idea of automatic differentiation to higher derivatives?
Let's check if we can repeat above Halley's method rate iterations.

Since we need 3 numbers, lets call it T (for triple) number.

Code:

T = {} -- note: f1, f2, f3 = f, f', f''

function T.new(f1,f2,f3) return setmetatable({f1, f2 or 0, f3 or 0}, T) end
function T.__add(f,g) return T.new(f[1]+g[1], f[2]+g[2], f[3]+g[3]) end
function T.__sub(f,g) return T.new(f[1]-g[1], f[2]-g[2], f[3]-g[3]) end
function T.__div(f,g) return T.rcp(g) * f end
function T.__unm(f) return T.new(-f[1], -f[2], -f[3]) end -- negate

-- (f*g)'  = f*g' + f'*g
-- (f*g)'' = f*g'' + f'*g' + f'*g' + f''*g

function T.__mul(f,g)
    local z = f[1]*g[3] + 2*f[2]*g[2] + f[3]*g[1]
    return T.new(f[1]*g[1], f[1]*g[2] + f[2]*g[1], z)
end

-- (1/f)'  = -f' / f^2
-- (1/f)'' = -(f^2 * f'' - (2*f*f') * f') / f^4

function T.rcp(f)
    return T.new(1/f[1], -f[2]/f[1]^2, (2*f[2]^2-f[1]*f[3]) / f[1]^3)
end

-- log(1+f)'  = f' / (1+f)
-- log(1+f)'' = ((1+f)*f'' - f'*f') / (1+f)^2

function T.log1p(f)
    local y = 1+f[1]
    return T.new(log1p(f[1]), f[2]/y, (y*f[3] - f[2]^2) / (y*y))
end

-- expm1(f)'  = exp(f) * f'
-- expm1(f)'' = exp(f) * f'' + exp(f)*f' * f'

function T.expm1(f)
    local x = expm1(f[1])
    local y, z = f[2]+f[2]*x, f[3]+f[3]*x
    return T.new(x, y, z + f[2]*y)
end

function T.halley(x,f)
    local r = -f[1]/f[2]
    r = r / (1 + 0.5*r*f[3]/f[2])
    x[1] = x[1] + r
    return x[1], f[1]
end

It work !

lua> n, pv, pmt, fv = T.new(10), T.new(50), T.new(-30), T.new(100)
lua> function f(x) return ((pv+fv)/T.expm1(T.log1p(x)*n) + pv)*x + pmt end
lua> x = T.new(pmt[1]/fv[1], 1) -- small edge
lua> for k=1,3 do print(T.halley(x, f(x))) end
-0.2844446981069045      1.3080888941077191
-0.28443599888025756    0.0007214129592014729
-0.28443599888025595    1.3500311979441904e-013

Note: outputs are (extrapolated x, f(x)); Shift up 2nd column to get points (x, f(x))