Finding the truncation error in an infinite sequence

Hi there,

Truncation error happens when an infinite series is ignored except for a small subset of its values. For instance, the exponential function e^x could be written as the infinite series total of 1 + x + x^2/2 + x^3/6 +... + x^n/n! +....

Any finite number of n will result in an approximation of the value of e^x that is inaccurate, but by increasing n, the error can be reduced to the desired level but this increase computation time.

How can I speed up the calculation by determining the dynamic level where error can be reduced to the barest minimum?

Any assistance would be greatly appreciated.

Many thanks

[James Tursa]

Do you have a specific series in mind?

[Life is Wonderful]

You could consider a 6 order Taylor series.

[James Tursa]

There are an infinite number of Taylor series. I meant what is the specific series? E.g., for exp(x) or sin(x) or ...?

And over what domain of x? I.e., what specific problem are you working on?

If you are just asking generally then we probably can't give you a definitive answer. In particular, getting a relatively fast and accurate calculation of a function value might involve complicated domain reduction techniques and using rational functions instead of Taylor series.

[Life is Wonderful]

Thank you @James Tursa

I'm looking for exp(x) and would like to know which technologies you believe are making the exponential fast and accurate; if you could elaborate, I'd be grateful.

[Walter Roberson]

The following discusses range reduction of the base for calculating exponentiation.

https://stackoverflow.com/questions/32409022/how-can-i-calculate-exponential-using-cordic-for-numbers-outside-1-1

The context for why people were looking for that has to do with an algorithm called CORDIC. MATLAB already has a fixed-point CORDIC exponential.

https://www.mathworks.com/help/fixedpoint/ref/cordiccexp.html

To put it another way: people do not use taylor series for exp() for FPGA purposes, they use a different successive-approximation approach that is well suited for FPGA.

[Life is Wonderful]

Thank you @Walter Roberson,

insight and much-needed input. To me, it sounds good.

Answer 1

[Torsten]

How can I speed up the calculation by determining the dynamic level where error can be reduced to the barest minimum?

You mean you want to stop adding x^n / n! when the sum of the rest series is smaller than a prescribed value ? And you ask how you can do this ?

This can only be done theoretically beforehand without the use of MATLAB.

[Life is Wonderful]

@Torsten

I'm talking about the non-deterministic error threshold. A quick example is provided below.

format long e
[exp(log(2)/2), 1+log(2)/2 + (log(2)/2)^2/factorial(2) + (log(2)/2)^3/factorial(3) + (log(2)/2)^4/factorial(4),(exp(log(2)/2) - 1+log(2)/2 + (log(2)/2)^2/factorial(2) + (log(2)/2)^3/factorial(3) + (log(2)/2)^4/factorial(4))]

ans = 1×3

1.414213562373095e+00 1.414169363672077e+00 8.283829260451723e-01

[exp(-log(2)/2), 1+(-log(2)/2) + (-log(2)/2)^2/factorial(2) + (-log(2)/2)^3/factorial(3) + (-log(2)/2)^4/factorial(4),(exp(log(2)/2) - 1+log(2)/2 + (log(2)/2)^2/factorial(2) + (log(2)/2)^3/factorial(3) + (log(2)/2)^4/factorial(4))]

ans = 1×3

7.071067811865475e-01 7.071461559459267e-01 8.283829260451723e-01

The error is different for each input; if the input number is small, the error condition can be reached quickly, say 1e-6; if the input number is large, we need more terms in the polynomial equation, i.e. more iteration; how can I make the threshold dynamic so that I don't have a fixed error tolerance number?

What I want to know is how to make the threshold dynamic so that I can decrease the polynomial's order and lower the number of terms, which will speed up execution time.

[Torsten]

As I said: Estimate the magnitude of the remaining sum theoretically. From this, you can determine N such that the error is less than a given threshold if you sum up at least until N or further.

The value you get for N differs depending on the infinite series and the x for which you want to evaluate the series. If you are only interested in the exponential function and its series, let me know. For this series, there are simple estimates of sum_{k=N}^{Inf} x^k / k! .

[Life is Wonderful]

@Torsten,

I'm not interested in the exponential function or its series, no. I'm seeking for an optimisation technique to lower the number of iterations after the acceptable level of accuracy is achieved.

[Torsten]

But for a general series, how should it be possible to know when an acceptable level of accuracy is achieved ? Say you have a series that is 1/n^2 up to the 10000th term and the 10001st term is 1e4. How could you deduce a stopping criterion for this series ?

Again: The series must be known in advance and theoretically analyzed to deduce rules on how far to sum to keep a prescribed error tolerance. You can works with many empirical rules, e.g. stop summing when the absolute value of the n-th term of the series is less than a certain tolerance multiplied by the absolute value of the sum so far. But this can fail due to unpredictable behaviour of the series.

[Life is Wonderful]

Yes, it is the problem I am trying to overcome. There are only a few methods accessible, including central, forward, and reverse substitution. Up until today, I hadn't noticed the true benefit. therefore I was looking for help.

Currently, I am terminating the infinite series when the absolute truncation error is smaller than a specified tolerance of fixed value, let's say 1e-8.

[Torsten]

There are only a few methods accessible, including central, forward, and reverse substitution.

Are these methods to speed up the convergence of the series ?

Currently, I am terminating the infinite series when the absolute truncation error is smaller than a specified tolerance of fixed value, let's say 1e-8.

But how do you know the absolute truncation error unless you know the function the series represents or you analyzed the series theoretically in advance ?

[Life is Wonderful]

There are only a few methods accessible, including central, forward, and reverse substitution.

Are these methods to speed up the convergence of the series ?

When considering a large number, which should be constrained before moving on to iterative series, I am not noticing any appreciable improvement.

But how do you know the absolute truncation error unless you know the function the series represents or you analyzed the series theoretically in advance ?

There isn't much job in front of me right now; I still have things to do. On it, I'm working. If I have any work to share, I jot it down.

[Torsten]

Don't spend too much time on it. There is no solution for this problem unless you know the series and you can estimate the truncation error in advance.

[Walter Roberson]

Consider for example the harmonic series, sum of 1/n . The first term is 1 and you know that by 10^16 that subsequent terms are each going to be be less than 1e-16 and when added to the initial 1 in double precision mathematics will not change the result. Therefore if you look at the relative values of the existing sum and see if the proposed new term is less than eps() of the partial sum, you can predict that there would be no point going beyond 1e16 terms at worst, and you would predict that the true value of the infinite harmonic sum is the finite value that you had arrived at.

You would, however, be wrong. The infinite harmonic sum is infinite.

The harmonic sum is the discrete approximation of the integral of 1/x, and that integral to X is ln(X). But ln(infinity) is infinity, giving evidence that the discrete approximation to it is also infinite.

The sum of the first N terms of the harmonic series is finite, but the infinite sum is infinite.

Therefore unless you do a mathematical analysis of the function you cannot use the fact that terms have become small to determine whether you have converged within floating point error.

[Torsten]

I also thought of this example, but I think it was assumed in advance that the series converges.

[Walter Roberson]

Consider 1/(n*exp(-c*n)) where c is positive and small. This decays exponentially and so should converge. But by making c small enough like 1e-50 you can make the decay be smaller than eps for n as small as 1e16 so you should be able to get into a situation where enough terms individually smaller than eps(1) add up to a large value and yet still be convergent

[Walter Roberson]

Hmmm, looks like that doesn't converge either.

Still, I imagine there must be some function that would have those properties, of having a large trailing residue of small values.

[Torsten]

sum_{n=1}^{Inf} 1/(n^(1+eps))

for eps > 0 arbitrarily small, e.g.

For c < 0, your example also works because

sum_{n=1}^{Inf} 1/n * exp(-c*n) <= sum_{n=0}^{Inf} exp(-c*n) = 1/(1-exp(-c)) < Inf

for c > 0.

[Life is Wonderful]

@Torsten,

Sorry for the delayed response; I was not at work.

I had trouble comprehending; if you could please illustrate your equation with numbers so that I can grasp it, that would be helpful.

[Walter Roberson]

@Life is Wonderful

My claim is that there exists convergent formulas for which the individual terms go to less than eps() by some finite boundary, but that there is still a lot of value in the tail to infinity

for example, if the terms are 1/(n^(1+1e-10)) then the sum does converge, to zeta(10000000001/10000000000) which is about 9999999173.17357 but the sum to 1e16 should be only slightly different than eulergamma + psi(1e16 + 1) which is about 37.4 .

Consider that for each positive n, 1/(n^(1+1e-10)) must be strictly less than 1/n then it follows that the sum of the modified term to 1e16 must be strictly less than the sum of 1/n to 1e16, which is the 37.4-ish number. So the long tail past 1e16 must overwhelm the terms to 1e16

syms n
expr = 1/(n^(1+sym(1e-10)))

expr = 

symsum(expr, n, 1, inf)

ans = 

vpa(ans)

ans = 

10000000000.57721566490881444517

symsum(1/n, n, 1, 1e16)

ans = 

vpa(ans)

ans = 

37.418577152806263854894375365032

From all of this, it follows that if you do not have a mathematical analysis of the function in order to be able to determine the appropriate stopping point, then stopping at any particular term on the basis of its absolute value being "small" is not sufficient to get you to within a given floating point truncation error.

Your task is impossible, at least in the form you stated.

Now, there are certainly a large number of functions that you could meaningfully stop when the individual terms got "small enough", but unless you can restrict the permitted functions to process, then you cannot handle the general situation.

[Life is Wonderful]

Thank you very much.

Please offer your thoughts once more.

What are your thoughts on the following code, and do you believe the best deterministic error tolerance is possible? What method can assist me in determining order 8 for all possible values ranging from -3.465735902799726e-01 to 3.465735902799726e-01 ?

If you look closely at the analysis, you'll note that I don't require new terms beyond the sixth order, therefore I'd like a deterministic technique for determining when to stop adding more terms or when to add more terms to the series, or there's a chance that I don't need 6th order polynomial terms.

clear all;clc;
syms x
f = exp(x);
T6 = taylor(f,x);
T8 = taylor(f,x,'Order',8);
T10 = taylor(f,x,'Order',10);

x = [-3.465735902799726e-01,3.465735902799726e-01];
subs(T6);
subs(T8);
subs(T10);

% plot exponential orders
figure(2);
fplot([T6 T8 T10 f])
grid on
legend('approximation of exp(x) with error O(x^6)', ...
       'approximation of exp(x) with error O(x^8)', ...
       'approximation of exp(x) with error O(x^{10})', ...
       'exp(x)','Location','Best')
title('Taylor Series Expansion')

clear all;clc;
format long e
syms x
f = exp(x);
T6 = taylor(f,x)

T6 = 

x = [-3.465735902799726e-01,3.465735902799726e-01];
subs(T6);

[exp(-3.465735902799726e-01), exp(3.465735902799726e-01)]

ans = 1×2

7.071067811865476e-01 1.414213562373095e+00

vpa(ans)

ans = 

syms x
f = exp(x);
T8 = taylor(f,x, 'Order',8)

T8 = 

x = [-3.465735902799726e-01,3.465735902799726e-01];
subs(T8);

[exp(-3.465735902799726e-01), exp(3.465735902799726e-01)]

ans = 1×2

7.071067811865476e-01 1.414213562373095e+00

vpa(ans)

ans = 

syms x
f = exp(x);
T4 = taylor(f,x, 'Order',4)

T4 = 

x = [-3.465735902799726e-01, 3.465735902799726e-01];
subs(T4);

[exp(-3.465735902799726e-01), exp(3.465735902799726e-01)]

ans = 1×2

7.071067811865476e-01 1.414213562373095e+00

vpa(ans)

ans = 

[Torsten]

If you look closely at the analysis, you'll notice that I don't need new terms beyond the sixth order, thus I'd like a deterministic technique to determining when to stop adding more terms or add more terms to the series

You don't know the value of the function the series represents (in this case [exp(-3.465735902799726e-01), exp(3.465735902799726e-01)]). So you don't know when to stop adding.

[Life is Wonderful]

Yes, I don't know how to stop, add, or consider the number of terms equation because the input could be a random number, and there are multiple input values to consider.

Note 3.465735902799726e-01 I gave an example. This cannot be true in real time, so I must use a random number as an input for the exponential function. So a method is required to approximate which can be utilised in a deterministic way for the number of terms to be used, so I am describing how to select a dynamic tolerence value in an infinite series.

[Torsten]

Yes, I don't know how to stop, add, or consider the number of terms equation because the input could be a random number, and there are multiple input values to consider.

No, the x values are not the problem. But you don't know that your series converges to exp(x). Thus you don't know which function f you need to insert when you consider f(x) - T6(x) or f(x) - T8(x) to estimate the error you make.

[Life is Wonderful]

Hmm, I'm fine with having a precise tolerance of 1e-6.

f(x) - T6(x) or f(x) - T8(x), this is a centre and forward method, and I'm arguing that if the requisite precision is achieved, why would I do extra iterations, which cost FPGA more resources?

Alternatively, if you could assist me in completing a bitwise precision check, I would be grateful.

[Torsten]

I give up because you don't seem to understand the point.

You don't know and thus cannot evaluate f(x) to estimate the error you make.

[Life is Wonderful]

That's all right (f(x) - I'm not aware. So another option I'm considering is what method I can use to identify a closer solution, i.e. how to verify precise matches with a certain tolerance?

how to perform a precision comparison here

[Walter Roberson]

FPGA cannot work with function handles, so if you are going for FPGA implementation you either need to choose a single fixed function or else you need to find a way to encode functions as data.

Also, for FPGA implementation you would more commonly be switching to fixed point calculations using CORDIC algorithms, as floating point cores take a lot of space in FPGA.

[Life is Wonderful]

FPGA implementation you would more commonly be switching to fixed point calculations using CORDIC algorithms

Yes, I'll use fixed point for the final implementation. But before hand , I was asking and making sure I 've a dead end and there's isn't a solution that exsist.

I was reading about vpasolve and running several tests concurrently.

As a substitute, I'm considering the precision check up to a specific decimal point. Please provide any useful ideas you may have.

Many thanks

[Walter Roberson]

vpasolve cannot be compiled or deployed to hardware.

I am not sure: are you trying to create an FPGA implementation of exp(x) specifically, or are you trying to get this to work for multiple functions? If you are trying to get this to work for multiple functions, are you looking for a code flow that will accept a function and automatically generate code suitable for evaluating the function to given precision over a given (input) range of values?

[Life is Wonderful]

employing vpasolve I'm attempting to comprehend how it functions in various situations so that I can learn more.

Anyway, Since I am unaware of multiple functions, it would be great if you could assist me with a code flow that accepts and generates code.

Thank you

[Walter Roberson]

Consider the simple looking

symsum(1/n, n, 1, N) - 25

and find the smallest N that makes the value >= 0

Now try again but with a target of (25 - 1e-6)

Is it practical to do even the smaller of the number of iterations?

If you were to taylor the sum then you would end up with coefficients at least as large as factorial(4e11) which is a value that exceeds 10^(4e12).

[Walter Roberson]

Therefore with even relatively simple easy to understand functions with modest targets, there is no practical way to do what you want. Taylor series are not powerful enough to be useful to calculate many of the transcendental functions.

Try even just sin(x) and a need to cover more than one period. You will not be able to get it to work.

[Life is Wonderful]

Is this the correct way to proceed?

syms n N
F(n,N) = symsum(1/n, n, 1, N) - 25;
F(n,N) = [3.465735902799726e-01,8];
subs(F(n,N))

ans = 

vpa(ans)

ans = 

syms n N
F(n,N) = symsum(1/n, n, 1, N) - (25  - 1e-6);
subs(F(n,N))

ans = 

vpa(ans)

ans = 

syms n N
F(n,N) = symsum(1/n, n, 1, N) - (25 - 1e-6);
F(n,N) = [0,10];
subs(F(n,N))

ans = 

[Walter Roberson]

syms n N
eqn1 = symsum(1/n, n, 1, N) - 25

eqn1 = 

N1 = ceil(vpasolve(eqn1))

N1 = 

40427833596

eqn2 = symsum(1/n, n, 1, N) - (25 - 1e-6)

eqn2 = 

N2 = ceil(vpasolve(eqn2))

N2 = 

40427793168

N1 - N2

ans = 

40428

So you would have to sum 40427833596 terms of 1/n to reach at least 25, and the number needed to reach within 1e-6 of 25 is 40427793168, which is about 40000 fewer.

Realistically speaking, you are never going to count up to even 40427793168 .

At what point do the individual terms become less than eps? That would obviously be at

uint64(1/eps())

ans = uint64 4503599627370496

which is even less likely for you to count to.

Since your tolerance is 1e-6, can you just stop at the point the terms become less than 1e-6 ?

vpa(subs(eqn1, N, 1e6+1))

ans = 

No, stopping at 1e6 iterations gives you a result that is only about 60% of the needed value.

[Life is Wonderful]

Thank you @Walter Roberson

I'm curious where the number 25 originates from ?

[Walter Roberson]

The 25 was arbitrary

[Life is Wonderful]

Thank you, the idea is now clear to me. Will C/C++ and HDL code be supported by the codegeneration process, which is the next step?

[Walter Roberson]

Nothing in the symbolic toolbox can be compiled. If you were to use the symbolic taylor series then you would need to matlabFunction() in a different session asking to write to file, and then in the other section to be compiled you would invoke the function by file name.

Note that if you go beyond 15 terms of taylor then the factorial(15) will exceed double precision resolution and your calculation will become unreliable. This gets to be a serious problem for transcendental functions.

[Life is Wonderful]

Yes, I am aware that factorial 20 would be the final one since, if your data type has a 64-bit width, it cannot fit beyond this amount.

But I'm not doing anything on the implementation front.