Path: csiph.com!x330-a1.tempe.blueboxinc.net!usenet.pasdenom.info!aioe.org!news.stack.nl!newsfeed.xs4all.nl!newsfeed5.news.xs4all.nl!xs4all!post.news.xs4all.nl!not-for-mail
Return-Path: <python-python-list@m.gmane.org>
X-Original-To: python-list@python.org
Delivered-To: python-list@mail.python.org
X-Spam-Status: OK 0.000
X-Spam-Evidence: '*H*': 1.00; '*S*': 0.00; 'else:': 0.03; 'subject:Python': 0.05; '"""': 0.07; 'calculating': 0.07; 'function,': 0.07; 'something,': 0.07; 'python': 0.08; 'assert': 0.09; 'delta': 0.09; 'function:': 0.09; 'mentions': 0.09; 'received:80.91': 0.09; 'received:80.91.229': 0.09; 'received:gmane.org': 0.09; 'received:list': 0.09; 'am,': 0.12; 'def': 0.13; 'operations,': 0.15; '"""use': 0.16; 'bouncing': 0.16; 'enough?': 0.16; 'err': 0.16; 'intrinsic': 0.16; 'otoh,': 0.16; 'skip:8 20': 0.16; 'thereto': 0.16; 'values:': 0.16; 'wrote:': 0.18; 'this?': 0.19; "haven't": 0.20; '(most': 0.21; 'trying': 0.21; 'discussion': 0.22; 'header:In-Reply-To:1': 0.22; 'feb': 0.22; 'numeric': 0.23; 'math': 0.24; 'traceback': 0.24; 'performing': 0.26; 'function': 0.27; 'import': 0.27; 'looks': 0.27; 'raise': 0.28; 'bit': 0.28; "i'm": 0.28; 'example': 0.29; 'problem': 0.29; 'print': 0.29; 'correct': 0.29; 'fairly': 0.30; '(the': 0.30; 'error': 0.30; 'anyone': 0.32; 'does': 0.32; 'break': 0.32; 'initial': 0.32; 'there': 0.33; 'this.': 0.33; 'header:User-Agent:1': 0.33; 'it.': 0.33; 'file': 0.34; 'header:X -Complaints-To:1': 0.34; 'from:addr:yahoo.co.uk': 0.34; 'mostly': 0.34; 'steven': 0.34; 'rather': 0.34; 'last):': 0.34; 'routine': 0.34; 'skip:i 40': 0.34; 'stuck': 0.34; 'try:': 0.34; 'operations': 0.35; 'to:addr:python-list': 0.35; 'none': 0.35; 'something': 0.35; 'be.': 0.35; 'apply': 0.35; 'two': 0.36; 'received:org': 0.36; 'with.': 0.37; 'but': 0.37; 'sometimes': 0.38; 'skip:d 40': 0.39; 'except': 0.39; 'point': 0.40; 'to:addr:python.org': 0.40; 'back': 0.60; 'your': 0.61; 'subject': 0.61; 'improve': 0.62; 'kind': 0.62; '19,': 0.68; 'received:89': 0.69; 'skip:2 20': 0.71; 'papers': 0.74; 'forth': 0.77; 'iterative': 0.84; 'one).': 0.84; 'received:as13285.net': 0.84; 'tricky': 0.84; 'suspected': 0.93; 'technique': 0.93
X-Injected-Via-Gmane: http://gmane.org/
To: python-list@python.org
From: Mark Lawrence <breamoreboy@yahoo.co.uk>
Subject: Re: Numeric root-finding in Python
Date: Sun, 12 Feb 2012 13:39:38 +0000
References: <4f375f0f$0$29986$c3e8da3$5496439d@news.astraweb.com> <e012b2f1-5754-443f-8a3f-a61528b3b8bc@m5g2000yqk.googlegroups.com>
Mime-Version: 1.0
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Content-Transfer-Encoding: 7bit
X-Gmane-NNTP-Posting-Host: host-89-243-199-55.as13285.net
User-Agent: Mozilla/5.0 (Windows NT 6.0; rv:10.0) Gecko/20120129 Thunderbird/10.0
In-Reply-To: <e012b2f1-5754-443f-8a3f-a61528b3b8bc@m5g2000yqk.googlegroups.com>
X-Antivirus: avast! (VPS 120212-0, 12/02/2012), Outbound message
X-Antivirus-Status: Clean
X-BeenThere: python-list@python.org
X-Mailman-Version: 2.1.12
Precedence: list
List-Id: General discussion list for the Python programming language <python-list.python.org>
List-Unsubscribe: <http://mail.python.org/mailman/options/python-list>, <mailto:python-list-request@python.org?subject=unsubscribe>
List-Archive: <http://mail.python.org/pipermail/python-list>
List-Post: <mailto:python-list@python.org>
List-Help: <mailto:python-list-request@python.org?subject=help>
List-Subscribe: <http://mail.python.org/mailman/listinfo/python-list>, <mailto:python-list-request@python.org?subject=subscribe>
Newsgroups: comp.lang.python
Message-ID: <mailman.5724.1329053992.27778.python-list@python.org>
Lines: 135
NNTP-Posting-Host: 2001:888:2000:d::a6
X-Trace: 1329053992 news.xs4all.nl 6962 [2001:888:2000:d::a6]:49505
X-Complaints-To: abuse@xs4all.nl
Xref: x330-a1.tempe.blueboxinc.net comp.lang.python:20270

On 12/02/2012 10:10, Eelco wrote:
> On Feb 12, 7:41 am, Steven D'Aprano<steve
> +comp.lang.pyt...@pearwood.info>  wrote:
>> This is only peripherally a Python problem, but in case anyone has any
>> good ideas I'm going to ask it.
>>
>> I have a routine to calculate an approximation of Lambert's W function,
>> and then apply a root-finding technique to improve the approximation.
>> This mostly works well, but sometimes the root-finder gets stuck in a
>> cycle.
>>
>> Here's my function:
>>
>> import math
>> def improve(x, w, exp=math.exp):
>>      """Use Halley's method to improve an estimate of W(x) given
>>      an initial estimate w.
>>      """
>>      try:
>>          for i in range(36):  # Max number of iterations.
>>              ew = exp(w)
>>              a = w*ew - x
>>              b = ew*(w + 1)
>>              err = -a/b  # Estimate of the error in the current w.
>>              if abs(err)<= 1e-16:
>>                  break
>>              print '%d: w= %r err= %r' % (i, w, err)
>>              # Make a better estimate.
>>              c = (w + 2)*a/(2*w + 2)
>>              delta = a/(b - c)
>>              w -= delta
>>          else:
>>              raise RuntimeError('calculation failed to converge', err)
>>      except ZeroDivisionError:
>>          assert w == -1
>>      return w
>>
>> Here's an example where improve() converges very quickly:
>>
>> py>  improve(-0.36, -1.222769842388856)
>> 0: w= -1.222769842388856 err= -2.9158979924038895e-07
>> 1: w= -1.2227701339785069 err= 8.4638038491998997e-16
>> -1.222770133978506
>>
>> That's what I expect: convergence in only a few iterations.
>>
>> Here's an example where it gets stuck in a cycle, bouncing back and forth
>> between two values:
>>
>> py>  improve(-0.36787344117144249, -1.0057222396915309)
>> 0: w= -1.0057222396915309 err= 2.6521238905750239e-14
>> 1: w= -1.0057222396915044 err= -2.6521238905872001e-14
>> 2: w= -1.0057222396915309 err= 2.6521238905750239e-14
>> 3: w= -1.0057222396915044 err= -2.6521238905872001e-14
>> 4: w= -1.0057222396915309 err= 2.6521238905750239e-14
>> 5: w= -1.0057222396915044 err= -2.6521238905872001e-14
>> [...]
>> 32: w= -1.0057222396915309 err= 2.6521238905750239e-14
>> 33: w= -1.0057222396915044 err= -2.6521238905872001e-14
>> 34: w= -1.0057222396915309 err= 2.6521238905750239e-14
>> 35: w= -1.0057222396915044 err= -2.6521238905872001e-14
>> Traceback (most recent call last):
>>    File "<stdin>", line 1, in<module>
>>    File "<stdin>", line 19, in improve
>> RuntimeError: ('calculation failed to converge', -2.6521238905872001e-14)
>>
>> (The correct value for w is approximately -1.00572223991.)
>>
>> I know that Newton's method is subject to cycles, but I haven't found any
>> discussion about Halley's method and cycles, nor do I know what the best
>> approach for breaking them would be. None of the papers on calculating
>> the Lambert W function that I have found mentions this.
>>
>> Does anyone have any advice for solving this?
>>
>> --
>> Steven
>
> Looks like floating point issues to me, rather than something
> intrinsic to the iterative algorithm. Surely there is not complex
> chaotic behavior to be found in this fairly smooth function in a +/-
> 1e-14 window. Otoh, there is a lot of floating point significant bit
> loss issues to be suspected in the kind of operations you are
> performing (exp(x) + something, always a tricky one).
>
> I would start by asking: How accurate is good enough? If its not good
> enough, play around the the ordering of your operations, try solving a
> transformed problem less sensitive to loss of significance; and begin
> by trying different numeric types to see if the problem is sensitive
> thereto to begin with.

HTH.

c:\Users\Mark\Python>type sda.py
import decimal
def improve(x, w, exp=decimal.Decimal.exp):
     """Use Halley's method to improve an estimate of W(x) given
     an initial estimate w.
     """
     try:
         for i in range(36):  # Max number of iterations.
             ew = exp(w)
             a = w*ew - x
             b = ew*(w + 1)
             err = -a/b  # Estimate of the error in the current w.
             if abs(err) <= 1e-16:
                 break
             print '%d: w= %r err= %r' % (i, w, err)
             # Make a better estimate.
             c = (w + 2)*a/(2*w + 2)
             delta = a/(b - c)
             w -= delta
         else:
             raise RuntimeError('calculation failed to converge', err)
         print '%d: w= %r err= %r' % (i, w, err)
     except ZeroDivisionError:
         assert w == -1
     return w
improve(decimal.Decimal('-0.36'), decimal.Decimal('-1.222769842388856'))
improve(decimal.Decimal('-0.36787344117144249'), 
decimal.Decimal('-1.0057222396915309'))
c:\Users\Mark\Python>sda.py
0: w= Decimal('-1.222769842388856') err= 
Decimal('-2.915897982757542086414504607E-7')
1: w= Decimal('-1.222770133978505953034526059') err= 
Decimal('-1.084120148360381932277303211E-19')
0: w= Decimal('-1.0057222396915309') err= 
Decimal('5.744538819905061986438230561E-15')
1: w= Decimal('-1.005722239691525155461180092') err= Decimal('-0E+2')

-- 
Cheers.

Mark Lawrence.