Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.python > #107579 > unrolled thread

Re: Optimizing Memory Allocation in a Simple, but Long Function

Back to article view | Back to comp.lang.python

This discussion starts older than the indexed window; earlier articles aren't shown. The article labeled Started by below is the oldest one visible, not the original post.

Contents

  Re: Optimizing Memory Allocation in a Simple, but Long Function Derek Klinge <schilke.60@gmail.com> - 2016-04-24 20:12 -0700
    Re: Optimizing Memory Allocation in a Simple, but Long Function Gregory Ewing <greg.ewing@canterbury.ac.nz> - 2016-04-25 19:39 +1200
      Re: Optimizing Memory Allocation in a Simple, but Long Function Oscar Benjamin <oscar.j.benjamin@gmail.com> - 2016-04-25 14:45 +0100
      Re: Optimizing Memory Allocation in a Simple, but Long Function Derek Klinge <schilke.60@gmail.com> - 2016-04-25 14:35 +0000
      Re: Optimizing Memory Allocation in a Simple, but Long Function Oscar Benjamin <oscar.j.benjamin@gmail.com> - 2016-04-26 10:54 +0100

#107579 — Re: Optimizing Memory Allocation in a Simple, but Long Function

Actually, I'm not trying to speed it up, just be able to handle a large
number of n.
(Thank you Chris for the suggestion to use xrange, I am on a Mac using the
stock Python 2.7)

I am looking at the number of iterations of linear approximation that are
required to get a more accurate representation.
My initial data suggest that to get 1 more digit of e (the difference
between the calculated and expected value falls under 10**-n), I need a
little more than 10 times the number of iterations of linear approximation.

I actually intend to compare these to other methods, including limit
definition that you provided, as well as the geometric series definition.

I am trying to provide some real world data for my students to prove the
point that although there are many ways to calculate a value, some are much
more efficient than others.

I tried your recommendation (Oscar) of trying a (1+1/n)**n approach, which
gave me very similar values, but when I took the difference between your
method and mine I consistently got differences of ~10**-15. Perhaps this is
due the binary representation of the decimals?

Also, it seems to me if the goal is to use the smallest value of n to get a
particular level of accuracy, changing your guess of N by doubling seems to
have a high chance of overshoot. I found that I was able to predict
relatively accurately a value of N for achieving a desired accuracy. By
this I mean, that I found that if I wanted my to be accurate to one
additional decimal place I had to multiply my value of N by approximately
10 (I found that the new N required was always < 10N +10).

Derek

On Sun, Apr 24, 2016 at 4:45 PM, Derek Klinge <schilke.60@gmail.com> wrote:

> Actually, I'm not trying to speed it up, just be able to handle a large
> number of n.
> (Thank you Chris for the suggestion to use xrange, I am on a Mac using the
> stock Python 2.7)
>
> I am looking at the number of iterations of linear approximation that are
> required to get a more accurate representation.
> My initial data suggest that to get 1 more digit of e (the difference
> between the calculated and expected value falls under 10**-n), I need a
> little more than 10 times the number of iterations of linear approximation.
>
> I actually intend to compare these to other methods, including limit
> definition that you provided, as well as the geometric series definition.
>
> I am trying to provide some real world data for my students to prove the
> point that although there are many ways to calculate a value, some are much
> more efficient than others.
>
> Derek
>
> On Sun, Apr 24, 2016 at 2:55 PM, Oscar Benjamin <
> oscar.j.benjamin@gmail.com> wrote:
>
>> 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
>> --
>> https://mail.python.org/mailman/listinfo/python-list
>>
>
>

[toc] | [next] | [standalone]

#107586

Derek Klinge wrote:
> Also, it seems to me if the goal is to use the smallest value of n to get a
> particular level of accuracy, changing your guess of N by doubling seems to
> have a high chance of overshoot.

If you want to find the exact n required, once you overshoot
you could use a binary search to narrow it down.

-- 
Greg

[toc] | [prev] | [next] | [standalone]

#107603

On 25 April 2016 at 08:39, Gregory Ewing <greg.ewing@canterbury.ac.nz> wrote:
> Derek Klinge wrote:
>>
>> Also, it seems to me if the goal is to use the smallest value of n to get
>> a
>> particular level of accuracy, changing your guess of N by doubling seems
>> to
>> have a high chance of overshoot.
>
>
> If you want to find the exact n required, once you overshoot
> you could use a binary search to narrow it down.

Also you can calculate the truncation error for Euler's method. Since

   f(t) = f(t0) + f'(t0)*(t - t0) + (1/2)f''(t0)*(t - t0)**2 + O((t - t0)**3)

Euler's method just uses the first two terms so

    x[n+1] = x[n] + dt*f(x[n])

the next term would be

   (1/2)*f'(x[n])*dt**2

Since in your case f'(x) = x and dt = 1/N that's

    (1/2)*x[n]*(1/N)**2

As a relative error (divide by x[n]) that's

    (1/2)*(1/N)**2

Let's add the relative error from N steps to get

    N*(1/2)*(1/N)**2 = 1/(2*N)

So the relative error integrating from 0 to 1 with N steps is 1/(2*N).
If we want a relative error of epsilon then the number of steps needed
is 1/(2*epsilon).

That is to say that for a relative error of 1e-4 we need N =
1/(2*1e-4) = 1e4/2 = 5e3 = 5000.

>>> import math
>>> N = 5000
>>> error = math.e - (1 + 1.0/N)**N
>>> relative_error = error / math.e
>>> relative_error
9.998167027596845e-05

Which is approximately 1e-4 as required.

--
Oscar

[toc] | [prev] | [next] | [standalone]

#107605

A couple thoughts. I think my original approach would be faster than binary
search for finding the minimum value of N needed to get a decimal level of
absolute accuracy from Euler's number. Here is my reasoning:
EulerlersNumber(13).LimitMethod() - math.e < .1 and
EulersNumber(135).LimitMethod - math.e < .01. N increased by a little more
than a factor of 10. My method would use 130, 131, 132, 133, 134, and 135
as guesses after 13. Using the double + binary search the guesses would be
26, 52, 104, 208, 156, 130, 143, 137, 134, 136, and then 136 after using
the information that n=13 gives an tolerance of < .1. If the trend holds,
then to get the next decimal of accuracy, you would need to do less than 10
searches each time. When I was using my linear approximation method I found
that I never had to do more than 10 additional guesses of N to get the next
digit of accuracy. When I was using EulersMethod() I found this trend held
for as long as I would allow my computer to do the calculations (though
accuracy to 10**-10).

What I find very interesting is that although the limit method is
mathematically equivalent to the linear approximation method, they give
different results, given sufficiently high values of N (at least in
Python). The limit method does not follow my predictions as accurately,
which leads to the question of whether or not the phenomenon I observed was
an artifact of rounding error or not.

Also, my explicit goal is to be able to handle large numbers of N, and to
reduce rounding error to a minimum. Using the fractions module to perform
the limit method with rational numbers rather than binary represented
decimals and got an error that the integer was too long to convert to float
and even when using smaller values of n (10**14) I seem to get values that
are not the same as when using decimals. It seems to me (I'm just thinking
about this at a low level and haven't written out any math to justify it
yet) that because the linear approximation method only multiplies and adds
the rounding error is smaller than in the limit method where numbers are
being taken to exponents that have rounding errors.

Although I see the value of relative error, I am just as interested in
absolute error (though admittedly they are directly related values).

Are there modules or libraries I can/should use to minimize rounding error
and use very large values of N and get an accurate answer? How does the
math module calculate the vale of e?

Thanks,
Derek

On Mon, Apr 25, 2016 at 6:49 AM Oscar Benjamin <oscar.j.benjamin@gmail.com>
wrote:

> On 25 April 2016 at 08:39, Gregory Ewing <greg.ewing@canterbury.ac.nz>
> wrote:
> > Derek Klinge wrote:
> >>
> >> Also, it seems to me if the goal is to use the smallest value of n to
> get
> >> a
> >> particular level of accuracy, changing your guess of N by doubling seems
> >> to
> >> have a high chance of overshoot.
> >
> >
> > If you want to find the exact n required, once you overshoot
> > you could use a binary search to narrow it down.
>
> Also you can calculate the truncation error for Euler's method. Since
>
>    f(t) = f(t0) + f'(t0)*(t - t0) + (1/2)f''(t0)*(t - t0)**2 + O((t -
> t0)**3)
>
> Euler's method just uses the first two terms so
>
>     x[n+1] = x[n] + dt*f(x[n])
>
> the next term would be
>
>    (1/2)*f'(x[n])*dt**2
>
> Since in your case f'(x) = x and dt = 1/N that's
>
>     (1/2)*x[n]*(1/N)**2
>
> As a relative error (divide by x[n]) that's
>
>     (1/2)*(1/N)**2
>
> Let's add the relative error from N steps to get
>
>     N*(1/2)*(1/N)**2 = 1/(2*N)
>
> So the relative error integrating from 0 to 1 with N steps is 1/(2*N).
> If we want a relative error of epsilon then the number of steps needed
> is 1/(2*epsilon).
>
> That is to say that for a relative error of 1e-4 we need N =
> 1/(2*1e-4) = 1e4/2 = 5e3 = 5000.
>
> >>> import math
> >>> N = 5000
> >>> error = math.e - (1 + 1.0/N)**N
> >>> relative_error = error / math.e
> >>> relative_error
> 9.998167027596845e-05
>
> Which is approximately 1e-4 as required.
>
> --
> Oscar
> --
> https://mail.python.org/mailman/listinfo/python-list
>

[toc] | [prev] | [next] | [standalone]

#107646

On 25 April 2016 at 15:35, Derek Klinge <schilke.60@gmail.com> wrote:
>
> Although I see the value of relative error, I am just as interested in
> absolute error (though admittedly they are directly related values).

I was referring to relative error because the relative error is the
same at each step making the calculation of the global error easier.

> Are there modules or libraries I can/should use to minimize rounding error
> and use very large values of N and get an accurate answer?

from decimal import Decimal, localcontext

def e(prec):
    with localcontext() as ctx:
        ctx.prec = 2*prec + 1
        N = Decimal(10)**(prec)
        eapprox = (1 + 1/N)**N
        ctx.prec = prec
        return +eapprox

print(e(50))

Alternatively you can just use Decimal(1).exp().

> How does the math
> module calculate the vale of e?

https://hg.python.org/cpython/file/tip/Include/pymath.h#l54

--
Oscar

[toc] | [prev] | [standalone]

Back to top | Article view | comp.lang.python

csiph-web