8
$\begingroup$

I am looking for a simple, inexpensive and very accurate approximation of the Lambert function ($W_0$ branch) ($-1/e < x < 0$).

$\endgroup$
7
$\begingroup$

Here is a post that I made on sci.math a while ago, regarding a method I have used.

Analysis of $we^w$

For $w\gt0$, $we^w$ increases monotonically from $0$ to $\infty$. When $w\lt0$, $we^w$ is negative.

Thus, for $x\gt0$, $\mathrm{W}(x)$ is positive and well-defined and increases monotonically.

For $w\lt0$, $we^w$ reaches a minimum of $-1/e$ at $w=-1$. On $(-1,0)$, $w e^w$ increases monotonically from $-1/e$ to $0$. On $(-\infty,-1)$, $w e^w$ decreases monotonically from $0$ to $-1/e$. Thus, on $(-1/e,0)$, $\mathrm{W}(x)$ can have one of two values, one in $(-1,0)$ and another in $(-\infty,-1)$. The value in $(-1,0)$ is called the principal value.

The iteration

Using Newton's method to solve $we^w$ yields the following iterative step for finding $\mathrm{W}(x)$: $$ w_{\text{new}}=\frac{xe^{-w}+w^2}{w+1} $$

Initial values of $w$

For the principal value, when $-1/e\le x\lt0$, and when $0\le x\le10$, use $w=0$. When $x\gt10$, use $w=\log(x)-\log(\log(x))$.

For the non-principal value, when $x$ is in $[-1/e,-.1]$, use $w=-2$; and if $x$ is in $(-.1,0)$, use $w=\log(-x)-\log(-\log(-x))$.

This says that for the branch you want, use the iteration with an initial value of $w=0$.

$\endgroup$
4
$\begingroup$

Corless, et al. 1996 is a good reference for implementing the Lambert $W$. One key to being efficient is a good initial guess for the iterative refinement. This is where series approximations are useful. For $(-1/\text{e}, 0)$ it's helpful to use two different expansions: about $0$ and about the branch point $-1/\text{e}$. Otherwise you can potentially waste a lot of time iterating.

Here's some unvectorized (for clarity) Matlab code for $W_0$. I've indicated the equation numbers from Corless, et al. on which various lines are based.

function w = fastW0(x)
% Lambert W function, upper branch, W_0
% Only valid for real-valued scalar x in [-1/e, 0]

if x > -exp(-sqrt(2))
    % Series about 0, Eq. 3.1
    w = ((((125*x-64)*x+36)*x/24-1)*x+1)*x;
else
    % Series about branch point, -1/e, Eq. 4.22
    p = sqrt(2*(exp(1)*x+1));
    w = p*(1-(1/3)*p*(1-(11/24)*p*(1-(86/165)*p*(1-(769/1376)*p))))-1;
end

tol = eps(class(x));
maxiter = 4;
for i = 1:maxiter
    ew = exp(w);
    we = w*ew;
    r = we-x; % Residual, Eq. 5.2
    if abs(r) < tol
        return;
    end

    % Halley's method, for Lambert W from Corless, et al. 1996, Eq. 5.9
    w = w-r/(we+ew-0.5*r*(w+2)/(w+1));
end

Here's a plot (for double precision, tolerance of machine epsilon) of number of the iterations until convergence if using no initial series guess, the series about $0$, the series about the branch point, or both series:

iterations vs. x

Of course the evaluation of the series approximation has a cost relative to performing an iteration. One can adjust how many terms are used in order to tune this.

$\endgroup$
3
$\begingroup$

An iterative procedure, given in wikipedia, q&d-translated to Pari/GP, which suits my needs well:

 LW(x, prec=1E-80, maxiters=200) = local(w, we, w1e); 
   w=0;    
   for(i=1,maxiters,  
           we=w*exp(w);
           w1e=(w+1)*exp(w); 
           if(prec>abs((x-we)/w1e),return(w));
           w = w-(we-x)/(w1e-(w+2)*(we-x)/(2*w+2)) ;
        );
    print("W doesn't converge fast enough for ",x,"( we=",we);
    return(0);

Example:

 default(precision,200)
   x = -0.99999/exp(1)
 %382 = -0.367875762377

   y=LW(x)
 %383 = -0.995534517079

  [chk=y*exp(y);x;chk-x]
 %384 = 
[-0.367875762377]
[-0.367875762377]
[2.49762470622 E-163]
$\endgroup$
  • $\begingroup$ The maximum number of iterations is : 1 ! $\endgroup$ – Claude Leibovici Jun 14 '13 at 11:00
  • $\begingroup$ @CLaude: ? I just tried x=-0.2 and got the result... $\endgroup$ – Gottfried Helms Jun 14 '13 at 11:17
  • $\begingroup$ Hi Gottfried ! Try at x=-0.999/e $\endgroup$ – Claude Leibovici Jun 14 '13 at 11:19
  • $\begingroup$ @Claude: perhaps you use too low precision? I'm using 200 digits per default, so I'm possibly overlooking that problem? I've appended one run for checking in my answer. $\endgroup$ – Gottfried Helms Jun 14 '13 at 11:32
  • $\begingroup$ The question arises: is this an approximation or a calculation of $LambertW$? Maple produces $evalf(LambertW(- \frac 1 {\exp(1)}),20)= 1.$. $\endgroup$ – user64494 Jun 14 '13 at 11:49
3
$\begingroup$

Maple produces $$\operatorname{Lambert W}(x)=x-{x}^{2}+{\frac {3}{2}}{x}^{3}-{\frac {8}{3}}{x}^{4}+{\frac {125}{24 }}{x}^{5}+O \left( {x}^{6} \right), x\to 0.$$ Let us compare $\operatorname{Lambert W}(-.2)=-.2591711018$ with $-0.2-(-0.2)^{2}+3/2(-0.2)^{3}-8/3(-0.2)^{4}+{\frac {125}{24}}(-0.2)^{5}= -.2579333334$.

$\endgroup$
  • $\begingroup$ Maple also produces its Pade approximation $$LambertW(x)= \left( {\frac {451}{340}}\,{x}^{3}+{\frac {228}{85}}\,{x}^{2}+x \right) \left( 1+{\frac {313}{85}}\,x+{\frac {1193}{340}}\,{x}^{2}+{ \frac {133}{204}}\,{x}^{3} \right) ^{-1} $$ which gives the value $LambertW(-0.2)=-.2591579777$. $\endgroup$ – user64494 Jun 14 '13 at 9:25
  • $\begingroup$ Thanks for this info. The problem is that I need the approximation to be almost as good close to zero as close to -1/e. $\endgroup$ – Claude Leibovici Jun 14 '13 at 10:58
  • $\begingroup$ The Pade approximation $[3,3]$ is that you want:see a pdf file. $\endgroup$ – user64494 Jun 14 '13 at 11:32
2
$\begingroup$

For the range $-\frac 1e \leq x \leq 0$, using $y=\sqrt{2(1+ex)}$ (as in this post) and building Padé approximants $P_n$ around $y=0$, what is obtained is $$P_1=\frac{-1+\frac{2 y}{3} } {1+\frac{y}{3} }$$ $$P_2=\frac{-1+\frac{14 y}{45}+\frac{301 y^2}{1080} } {1+\frac{31 y}{45}+\frac{83 y^2}{1080} }$$ $$P_3=\frac{-1-\frac{974 y}{22659}+\frac{21865 y^2}{51792}+\frac{89131 y^3}{1165320} } {1+\frac{23633 y}{22659}+\frac{104225 y^2}{362544}+\frac{3167 y^3}{196560} }$$ $$P_4=\frac{-1-\frac{11637254 y}{29330279}+\frac{463636649 y^2}{1055890044}+\frac{23930361857 y^3}{110868454620}+\frac{192684057311 y^4}{10643371643520} } {1+\frac{40967533 y}{29330279}+\frac{659231191 y^2}{1055890044}+\frac{1928737771 y^3}{20157900840}+\frac{34384971553 y^4}{10643371643520} }$$ The maximum errors occur for $x=0$; for the above formulae, they are $\Delta_1=-3.88\times 10^{-2}$, $\Delta_2=-1.23\times 10^{-3}$, $\Delta_3=-3.73\times 10^{-5}$ and $\Delta_4=-1.12\times 10^{-6}$.

Concerning $$S_i=\int_{-\frac 1e}^0 \left(W(x)-P_i\right)^2\,dx$$ they are $S_1=3.93\times 10^{-4}$, $S_2=2.57\times 10^{-7}$, $S_3=1.76\times 10^{-10}$ and $S_4=1.26\times 10^{-13}$.

Edit

For a "sanity" check, I also considered all the possible $[k,6-k]$ Padé approximants and obtained the following results $$\left( \begin{array}{ccc} k & \Delta_k & S_k \\ 0 & -0.108806 & 1.84 \times 10^{-3} \\ 1 & -0.007522 & 6.16 \times 10^{-6} \\ 2 & -0.000171 & 3.49 \times 10^{-9}\\ \color{red} { 3} & \color{red} {-0.000037} & \color{red} {1.76 \times 10^{-10}}\\ 4 & -0.000082 & 8.17 \times 10^{-10}\\ 5 & -0.000966 & 1.05 \times 10^{-7}\\ 6 & -0.094985 & 9.04 \times 10^{-4} \end{array} \right)$$

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.