R help - Nov 2003 - Lambert's W function

Hello List

does anyone have an R function for the Lambert W function?  I need 
complex arguments.

[the Lamert W function W(z) satisfies

W(z)*exp(W(z)) = z

but I could'nt even figure out how to use uniroot() for complex z]





-- 
Robin Hankin
Uncertainty Analyst
Southampton Oceanography Centre
SO14 3ZH
tel +44(0)23-8059-7743
initialDOTsurname at soc.soton.ac.uk (edit in obvious way; spam precaution)

On Tue, 2003-11-25 at 12:09, Robin Hankin wrote:
> Hello List
> 
> does anyone have an R function for the Lambert W function?  I need 
> complex arguments.
> 
> [the Lamert W function W(z) satisfies
> 
> W(z)*exp(W(z)) = z
> 
> but I could'nt even figure out how to use uniroot() for complex z]
> 
> 
There are several 'branches' of W(z) for complex arguments, since the
inverse of f(z) = z * exp(z) does not exist. Therefore you have to
constrain root-finding methods appropriately.

Here is an implementation for real valued arguments, if that is of any
help to you:

WE1 <- function(x) {
	if (1) {
	o <- optim(c(1),
		# optimize squared difference of z and log(t) + t
		function(t) { (x - (log(t[1]) + t[1]))**2 },
		# the gradient of the above function
		#function(t) { 2(log(t[1]) + t[1] - x)(1/t[1] - 1) },
		method="L-BFGS-B", lower=exp(-70),
			control = list(lmm=40, maxit=10000, factr=1e9));
	} else {
	o <- optim(c(1),
		# optimize squared difference of z and log(t) + t
		function(t) { (x - (log(t[1]) + t[1]))**2 },
		method="BFGS");
	}
	o$par
}

W2 <- function(x) {
	eps <- 10^-5;	# precision

	if (x == 0) {
		wnew <- 0
	} else {
		# initial guess
		wold <- if (x < 2.5) {
			(x + 4/3*x^2)/(1 + 7/3*x + 5/6*x^2)
		} else {
			log(x)
		};
		if (x < 0) print("error");
			lx <- log(x);
		wnew <- 2 * wold;
		if (x != Inf) {
			while (abs( (wnew - wold)/wold ) > eps) {
				wold <- wnew
				zn <- lx - wold - log(wold);
				y <- 2*(1 + wold)*((1 + wold) + 2/3*zn) - zn;
				wnew <- wold*(1 + (zn * y)/((1 + wold)*(y - zn)));
			}
		} else {
			wnew <- Inf;
		}
	}
	wnew
}
WE <- function(z) {
	if ( z < -1 ) { W2(exp(z)) }
	else { WE1(z) };
}
W <- function(z) { WE(log(z)) }


Stefan

This implementation, originally written by Nici Schraudolph, allows you
to choose which branch you want.  I've checked the answers for complex
arguments, non-systematically, against Mathematica's ProductLog function.

  Ben Bolker

lambertW = function(z,b=0,maxiter=10,eps=.Machine$double.eps,
  min.imag=1e-9) {
  if (any(round(Re(b)) != b))
    stop("branch number for W must be an integer")
  if (!is.complex(z) && any(z<0)) z=as.complex(z)
  ## series expansion about -1/e
  ##
  ## p = (1 - 2*abs(b)).*sqrt(2*e*z + 2);
  ## w = (11/72)*p;
  ## w = (w - 1/3).*p;
  ## w = (w + 1).*p - 1
  ##
  ## first-order version suffices:
  ##
  w = (1 - 2*abs(b))*sqrt(2*exp(1)*z + 2) - 1
  ## asymptotic expansion at 0 and Inf
  ##
  v = log(z + as.numeric(z==0 & b==0)) + 2*pi*b*1i;
  v = v - log(v + as.numeric(v==0))
  ## choose strategy for initial guess
  ##
  c = abs(z + exp(-1));
  c = (c > 1.45 - 1.1*abs(b));
  c = c | (b*Im(z) > 0) | (!Im(z) & (b == 1))
  w = (1 - c)*w + c*v
  ## Halley iteration
  ##
  for (n in 1:maxiter) {
    p = exp(w)
    t = w*p - z
    f = (w != -1)
    t = f*t/(p*(w + f) - 0.5*(w + 2.0)*t/(w + f))
    w = w - t
    if (abs(Re(t)) < (2.48*eps)*(1.0 + abs(Re(w)))
        && abs(Im(t)) < (2.48*eps)*(1.0 + abs(Im(w))))
      break
  }
  if (n==maxiter) warning(paste("iteration limit (",maxiter,
        ") reached, result of W may be inaccurate",sep=""))
  if (all(Im(w)<min.imag)) w = as.numeric(w)
  return(w)
}

---------

\name{lambertW}
\alias{lambertW}
\title{Lambert W function}
\description{
  Computes the Lambert W function, giving efficient solutions to the equation
x*exp(x)==x
}
\usage{
lambertW(z, b = 0, maxiter = 10, eps = .Machine$double.eps, min.imag = 1e-09)
}
\arguments{
  \item{z}{(complex) vector of values for which to compute the function}
  \item{b}{(integer) vector of branches: b=0 specifies the principal
    branch, 0 and -1 are the ones that can take non-complex values}
  \item{maxiter}{maximum numbers of iterations for convergence}
  \item{eps}{convergence tolerance}
  \item{min.imag}{maximum magnitude of imaginary part to chop when
    returning solutions}
}
\details{
Compute the Lambert W function of z.  This function satisfies
W(z)*exp(W(z)) = z, and can thus be used to express solutions
of transcendental equations involving exponentials or logarithms.
The Lambert W function is also available in 
Mathematica (as the ProductLog function), and in Maple.
}
\value{
  Complex or real vector of solutions.
}
\references{Corless, Gonnet, Hare, Jeffrey, and Knuth (1996), "On the
Lambert
W Function", Advances in Computational Mathematics 5(4):329-359}
\author{Nici Schraudolph <schraudo at inf.ethz.ch> (original
  version (c) 1998), Ben Bolker (R translation)
  }
\note{
This implementation should return values within 2.5*eps of its
counterpart in Maple V, release 3 or later.  Please report any
discrepancies to the author or translator.
}
\examples{
curve(lambertW(x),from=0,to=10)
pvec = seq(-1,1,length=40)
m = outer(pvec,pvec,function(x,y)Re(lambertW(x+y*1i)))
persp(pvec,pvec,m,
      theta=290,shade=0.5,zlab="lambertW")
num1 = uniroot(function(x) {x*exp(x)-1},lower=0,upper=1,tol=1e-9)
abs(lambertW(1)-num1$root)<1e-9
}
\keyword{math}

-- 
620B Bartram Hall                            bolker at zoo.ufl.edu
Zoology Department, University of Florida    http://www.zoo.ufl.edu/bolker
Box 118525                                   (ph)  352-392-5697
Gainesville, FL 32611-8525                   (fax) 352-392-3704