On Wednesday, January 22, 2014 3:30:53 AM UTC-8, Hans W Borchers wrote: > > Jason, very nice, and I will certainly look at DualNumbers.jl. But I have > doubts about your approach to the Lambert W function. >
Yes, my implementation was not meant for production--just to demonstrate that automatic differentiation is compatible with Newton's method. > See, lambert_newton(-1/e) returns -0.9999999434202022 and should return -1. > Below that, the real branch of Lambert W is not defined, but > lambert_newton(-0.5) returns 0.9706698296525587 though we are in the > imaginary branch. > > This is not so problematic as you could return NaN for those values. What > about parameters x > 20, for example > > julia> lambert_newton(25) > 3.1737928264162476 > > while the true value is 2.36015... And lambert_newton(1000) returns NaN, > true value 5.249603... That is why CASs such as Maple define the Lambert W > function through a Newton-Raphson procedure such as > > function lambert_W(x::Real) > if x < -1/exp(1); return NaN; end > w0 = 1.0 > w1 = w0 - (w0 * exp(w0) - x)/((w0 + 1) * exp(w0) - > (w0 + 2) * (w0 * exp(w0) - x)/(2 * w0 + 2)) > n = 1 > while abs(w1 - w0) > 1e-15 && n <= 20 > w0 = w1 > w1 = w0 - (w0 * exp(w0) - x)/((w0 + 1) * exp(w0) - > (w0 + 2) * (w0 * exp(w0) - x)/(2 * w0 + 2)) > n += 1 > end > return w1 > end > > Now the value of lambert_W(-1/e) = -0.9999999911436428 is slightly better > than before and for x >= 20 and much higher we get correct results. > You've used Halley's method here, which is an extension of Newton-Raphson that uses 2nd derivative information to get cubic convergence rather than the quadratic convergence. The one place that I wouldn't expect either of them to converge rapidly to high accuracy is -1/e because that is a branch point with an undefined derivative
