Path: csiph.com!fu-berlin.de!uni-berlin.de!not-for-mail
From: Oscar Benjamin <oscar.j.benjamin@gmail.com>
Newsgroups: comp.lang.python
Subject: Re: Optimizing Memory Allocation in a Simple, but Long Function
Date: Sun, 24 Apr 2016 22:55:34 +0100
Lines: 200
Message-ID: <mailman.59.1461534957.32212.python-list@python.org>
References: <CAGuRns8KhOFL-a1_a=0hMP=NbktoJw78ih3j4uUSbvkF8bf3pw@mail.gmail.com> <CAPTjJmrQfk20EuoZVMLqKv8F5r=YyKj0nr8iNPw7DWoCEYLmYA@mail.gmail.com> <CAGuRns86c4G5BHX_FF-OBnV_eRVKyngx_Dfp0u7hoHRT0WN++w@mail.gmail.com> <CAPTjJmq5g848-PMTgHOg7sawm+JgA5mSfwDHA-ZZgsmYr6yLCA@mail.gmail.com> <CAGuRns8LsThq3iV-bkdBbg9fZO7s61z4ggBy=YkGDqhaB_J4Xw@mail.gmail.com> <CAGuRns9uvHVBXP79X1Sp8vsW4JpWF-3g1mJoWHOq-cGog61cug@mail.gmail.com> <CAPTjJmpvtP5pj_NiJNXGLmu2QCqSYZNgZjKisRRMRTtzvGHKyQ@mail.gmail.com> <CAHVvXxSwb8=a+b5WQfd3KeXTe8VuL4RMnwu3TSNp8129FJHJGw@mail.gmail.com>
Mime-Version: 1.0
Content-Type: text/plain; charset=UTF-8
In-Reply-To: <CAPTjJmpvtP5pj_NiJNXGLmu2QCqSYZNgZjKisRRMRTtzvGHKyQ@mail.gmail.com>
Precedence: list
Xref: csiph.com comp.lang.python:107573

On 24 April 2016 at 19:21, Chris Angelico <rosuav@gmail.com> wrote:
> On Mon, Apr 25, 2016 at 4:03 AM, Derek Klinge <schilke.60@gmail.com> wrote:
>> Ok, from the gmail web client:
>
> Bouncing this back to the list, and removing quote markers for other
> people's copy/paste convenience.
>
> ## Write a method to approximate Euler's Number using Euler's Method
> import math
>
> class EulersNumber():
>     def __init__(self,n):
>         self.eulerSteps = n
>         self.e    = self.EulersMethod(self.eulerSteps)
>     def linearApproximation(self,x,h,d): # f(x+h)=f(x)+h*f'(x)
>         return x + h * d
>     def EulersMethod(self, numberOfSteps): # Repeate linear
> approximation over an even range
>         e = 1                                                    # e**0 = 1
>         for step in range(numberOfSteps):
>             e = self.linearApproximation(e,1.0/numberOfSteps,e) # if
> f(x)= e**x, f'(x)=f(x)
>         return e
>
>
> def EulerStepWithGuess(accuracy,guessForN):
>     n = guessForN
>     e = EulersNumber(n)
>     while abs(e.e - math.e) > abs(accuracy):
>         n +=1
>         e = EulersNumber(n)
>         print('n={} \te= {} \tdelta(e)={}'.format(n,e.e,abs(e.e - math.e)))
>     return e
>
>
> def EulersNumberToAccuracy(PowerOfTen):
>     x = 1
>     theGuess = 1
>     thisE = EulersNumber(1)
>     while x <= abs(PowerOfTen):
>         thisE = EulerStepWithGuess(10**(-1*x),theGuess)
>         theGuess = thisE.eulerSteps * 10
>         x += 1
>     return thisE
>
>
>> To see an example of my problem try something like EulersNumberToAccuracy(-10)

Now that I can finally see your code I can see what the problem is. So
essentially you want to calculate Euler's number in the following way:

e = exp(1) and exp(t) is the solution of the initial value problem
with ordinary differential equation dx/dt = x and initial condition
x(0)=1.

So you're using Euler's method to numerically solve the ODE from t=0
to t=1. Which gives you an estimate for x(1) = exp(1) = e.

Euler's method solves this by going in steps from t=0 to t=1 with some
step size e.g. dt = 0.1. You get a sequence of values x[n] where

   x[0] = x(0) = 1  # initial condition
   x[1] = x[0] + dt*f(x[0]) = x[0] + dt*x[0]
   x[2] = x[1] + dt*x[1] # etc.

In order to get to t=1 in N steps you set dt = 1/N. So simplifying
your code (all the classes and functions are just confusing the
situation here):

N = 1000
dt = 1.0 / N
x = 1
for n in range(N):
    x = x + dt*x
print(x)

When I run that I get:
2.71692393224

Okay that's great but actually you want to be able to set the accuracy
required and then steadily increase N until it's big enough to achieve
the expected accuracy so you do this:

import math

error = 1
accuracy = 1e-2

N = 1
while error > accuracy:
    dt = 1.0 / N
    x = 1
    for n in range(N):
        x = x + dt*x
    error = abs(math.e - x)
    N += 1
print(x)

But what happens here? You have a loop in a loop. The inner loop takes
n over N values. The outer loop takes N from 1 up to Nmin where Nmin
is the smallest value of N such that we achieve the desired accuracy.

This is a classic case of a quadratic performance algorithm. As you
make the accuracy smaller you're implicitly increasing Nmin. However
the algorithmic performance is quadratic in Nmin i.e. O(Nmin**2). The
problem is the nested loops. If you have an outer loop that increases
the length of an inner loop by 1 at each step then you have a
quadratic algorithm akin to:

# This loop is O(M**2)
for n in range(N):
    for N in range(M):
        # do stuff

To see that it is quadratic see:

    https://en.wikipedia.org/wiki/Triangular_number

The simplest fix here is to replace N+=1 with N*=2. Instead of
increasing the number of steps by one if the accuracy is not small
enough then you should double the number of steps. That will give you
an O(Nmin) algorithm.

    https://en.wikipedia.org/wiki/1/2_%2B_1/4_%2B_1/8_%2B_1/16_%2B_%E2%8B%AF

A better method is to do a bit of algebra before putting down the code:

    x[1] = x[0] + h*x[0] = x[0]*(1+h) = x[0]*(1+1/N) = (1+1/N)
    x[2] = x[1]*(1+1/N) = (1+1/N)**2
    ...
    x[n] = (1 + 1/n)**n

So doing the loop for Euler's method is equivalent to just writing:

    x = (1 + 1.0/N)**N

This considered as a sequence in N is well known as a sequence that
converges to e. In fact this is how the number e was first discovered:

    https://en.wikipedia.org/wiki/E_%28mathematical_constant%29#Compound_interest

Python can compute this much quicker than your previous version:

N = 1
for _ in range(40):
    N *= 2
    print((1 + 1.0/N) ** N)

Which runs instantly and gives:

2.25
2.44140625
2.56578451395
2.63792849737
2.67699012938
2.69734495257
2.70773901969
2.71299162425
2.71563200017
2.71695572947
2.71761848234
2.71795008119
2.71811593627
2.71819887772
2.71824035193
2.7182610899
2.71827145911
2.71827664377
2.71827923611
2.71828053228
2.71828118037
2.71828150441
2.71828166644
2.71828174745
2.71828178795
2.71828180821
2.71828181833
2.7182818234
2.71828182593
2.71828182719
2.71828182783
2.71828182814
2.7182818283
2.71828182838
2.71828182842
2.71828182844
2.71828182845
2.71828182845
2.71828182846
2.71828182846

So your method is computing the above numbers but in a slower way that
also has more potential for rounding error. The error here is 1e-13
for the last numbers in this sequence. But N=2**40 so your Euler
method would need approximately 10**12 iterations in your inner loop
to get the same result. That's going to be slow even if you don't use
a quadratic algorithm.

--
Oscar