I spent a little time looking into Padé approximants, including how well they work for root-finding. These are ratios of polynomials which are used as approximations to other functions. For instance, for both polynomials of order 2, the form is a + bx + cx2 ⁄ 1 + dx + ex2 . Note the 1 in the denominator is chosen to avoid the ambiguity of being able to multiply through by a constant.
Just as with the derivation of a Taylor series, one can determine the coefficients of the Padé approximate by equating the derivatives of the desired function to the derivatives of the approximate. For instance, a is always equal to F(0).
If the order of the denominator is 0 (i.e., the denominator is 1), then the Padé approximants are just the Taylor series for the function. I looked at the approximate of order (0,1) as well (a ⁄ 1 + dx, which is F2 ⁄ F - F′·x), but found that it converged worse that the Taylor series of order 1.
The useful one seems to be the approximant of order (1,1). Since it has three free coefficients, it is equivalent to using a Taylor series up to order 2 (the x2 term). In my experiments the Padé gave closer approximations than the Taylor series. It is also easier to work with (especially if using rational numbers) than the Taylor series since there is no squared term. For the same reason, it is easier to solve for x and use the approximation to solve for a root.
Newton's method for root-finding equates to writing the Taylor expansion for a function out to order 1, then inverting it to obtain an expression for x. Obviously using another term in the Taylor expansion would give more accuracy, but since that produces a quadratic equation, the solution for x has a square root and is very messy (and in many cases, it is the square root function itself we are trying to solve in the first place). Since the Padé approximate of order (1,1) only has linear terms of x, it can still be easily solved, and provides quicker convergence to the root than Newton's method.
Once I had written out the equation for the correction term using this method (which is (F′ ⁄ α + F″ ⁄ 2F′)-1, where α is the needed correction to the function value), I saw that it was the same as that given by Halley's method. What I find interesting is, if you read the Wikipedia article, you'll see the derivation of Halley's method is completely different than the process I went through. The convergence of the expression from two different approaches increases my confidence that it is a good method for root-finding. As an added bonus, when computing with rational terms, the size of the terms grows more slowly than from Newton's method.
Showing posts with label Math. Show all posts
Showing posts with label Math. Show all posts
Sunday, August 30, 2009
Sunday, August 16, 2009
Memoize
I recently came across the computer science concept of memoization. What's cool is that I had made some suggestions on designs at work that amount to the same thing.
In general, memoization refers to having a computer function remember the problems it's already solved, and just return the answer again if it recognizes the problem, rather than compute it again, as most software does. The common example is the factorial function. Since the factorial is the product of all the numbers up to the argument, on the way to, say 6!, you will compute all the factorials less than 6. Usually the function would loop over all the numbers and multiply them together. Instead, a memoized version would have stored all the results up to the highest argument it had previously seen, and then only compute those beyond that (storing the new results) if necessary.
One type of object-oriented class that I had proposed is a smart Angle class. Since angles can be represented in many different ways, I wanted to be able to get any version the calling software wanted, without knowing what version it had been stored in. In addition to radians and degrees, an angle can be (and often is in our work) stored as the combination of sine and cosine values. Besides being able to use an angle without knowing what format it was in originally, I suggested that the other representations should be left uncomputed, but once they were asked for, the result should be stored to avoid unnecessary further computations. Clearly this can been seen as memoization of the conversion.
In the general case, the "hit rate" of a memoized function may not be that high, depending on the number and valid range of the inputs. However, for an object whose data never changes once it is created, essentially one of the inputs has been fixed, raising the chance that the same result will be requested repeatedly. (In my design, each Angle object would be for a particular angle. If the value of the variable changed, a new Angle object would be created to store it.)
By the way, I'm fascinated with the idea of "smart numbers" in computer programs - numbers as objects (or other collections of data, if not object-oriented) which store more than just their value, or that store exact representations in alternate formats. The Angle example given, quadratic irrationals, storing numbers as rationals or continued fractions instead of decimal, and numbers that store their own error bounds, are all concepts that I've thought about over the years. If anyone knows where one of these concepts has been used in a real program, I'd be interested to hear about it.
In general, memoization refers to having a computer function remember the problems it's already solved, and just return the answer again if it recognizes the problem, rather than compute it again, as most software does. The common example is the factorial function. Since the factorial is the product of all the numbers up to the argument, on the way to, say 6!, you will compute all the factorials less than 6. Usually the function would loop over all the numbers and multiply them together. Instead, a memoized version would have stored all the results up to the highest argument it had previously seen, and then only compute those beyond that (storing the new results) if necessary.
One type of object-oriented class that I had proposed is a smart Angle class. Since angles can be represented in many different ways, I wanted to be able to get any version the calling software wanted, without knowing what version it had been stored in. In addition to radians and degrees, an angle can be (and often is in our work) stored as the combination of sine and cosine values. Besides being able to use an angle without knowing what format it was in originally, I suggested that the other representations should be left uncomputed, but once they were asked for, the result should be stored to avoid unnecessary further computations. Clearly this can been seen as memoization of the conversion.
In the general case, the "hit rate" of a memoized function may not be that high, depending on the number and valid range of the inputs. However, for an object whose data never changes once it is created, essentially one of the inputs has been fixed, raising the chance that the same result will be requested repeatedly. (In my design, each Angle object would be for a particular angle. If the value of the variable changed, a new Angle object would be created to store it.)
By the way, I'm fascinated with the idea of "smart numbers" in computer programs - numbers as objects (or other collections of data, if not object-oriented) which store more than just their value, or that store exact representations in alternate formats. The Angle example given, quadratic irrationals, storing numbers as rationals or continued fractions instead of decimal, and numbers that store their own error bounds, are all concepts that I've thought about over the years. If anyone knows where one of these concepts has been used in a real program, I'd be interested to hear about it.
Thursday, August 13, 2009
A Tale of Two Sequences
Last week I came home happy one day because of a simple little math expression. The Fibonacci sequence has the cool property that for each term, if a prior term's index is a divisor of the current index, then that term is also a divisor of the current term. So, for example, F15 is divisible by F3 and F5 (and an additional factor). More formally, a divides b implies Fa divides Fb. There are other sequences with this property, including 10n − 1 (or any bn − 1). But the Fibonacci sequence is a linear recurrence relation (each term is the sum of the two previous terms) and the sequence of one less than powers of ten is exponential. I figured there must be some way to characterize what types of sequences can have this property.
Well, I still don't know if there are other types of sequences with this property, but I do now know why these two types share it. I came across this relationship in a paper I found online, and as soon as I saw it I realized it could be generalized. They are actually the same type, and can each be expressed as the other. For instance, 10n − 1 can also be obtained using Fn = 11·Fn-1 − 10·Fn-2, and each linear recurrence sequence can also be made into a difference-of-powers sequence, though the base "b" may not be rational.
As with many techniques in math and physics, once you know there is an equivalence, it is trivial to set up a set of equations in that form and figure out the coefficients. So, one the one hand I feel that I should have tried it sooner, but on the other I feel that if it were widely known then I would have run across mention of it sooner. Regardless, I am happy to know a little bit more about a set of sequences which I find fascinating (more properties will get their own post one day).
Well, I still don't know if there are other types of sequences with this property, but I do now know why these two types share it. I came across this relationship in a paper I found online, and as soon as I saw it I realized it could be generalized. They are actually the same type, and can each be expressed as the other. For instance, 10n − 1 can also be obtained using Fn = 11·Fn-1 − 10·Fn-2, and each linear recurrence sequence can also be made into a difference-of-powers sequence, though the base "b" may not be rational.
As with many techniques in math and physics, once you know there is an equivalence, it is trivial to set up a set of equations in that form and figure out the coefficients. So, one the one hand I feel that I should have tried it sooner, but on the other I feel that if it were widely known then I would have run across mention of it sooner. Regardless, I am happy to know a little bit more about a set of sequences which I find fascinating (more properties will get their own post one day).
Wednesday, July 1, 2009
Rational Approximation to Roots Continued...
This morning I was able to complete the symbolic expression for the limit of the series. For large x, the expression (1 + a⁄(bx-c))x can be expanded using the binomial theorem to 1 + x ⋅ a⁄(bx-c) + 1⁄2 x (x-1) (a⁄(bx-c))2 + 1⁄6 x (x-1) (x-2) (a⁄(bx-c))3 + … ; only keeping the highest power of x in both the numerator and the denominator, they cancel and it becomes ∑ 1⁄n! (a⁄b)n.
This is where I was last night, but this morning I realized that it could be a MacLaurin series. A little browsing in the CRC book and sure enough, it is the series for ey, with y=a⁄b. So, the actual limit is ea⁄b, or solving for a⁄b, we want a⁄b = ln N; indeed, ln 2 is around .693, and 9/13 is a good approximation for it. Now that I know exactly what irrational number I am trying to have for the ratio a⁄b to converge to any given N, I can use continued fractions to choose good a and b to approximate that ratio.
By the way, to get an exact result at x=1, c should be set to b - a⁄N-1.
However, using Excel to plot the actual values as a function of x, it becomes clear that the convergence becomes slow for even N around 10. (7, 10, 3) gives values within .015 of 2, but the first few terms when converging to 10 are as low as 8! One can mitigate this somewhat by raising c such that the values for x=1 and x=2 are poor (do we really need an approximation to the 1st root?), bringing the errors for the next terms down, but this only reduces the error by about half.
So, if you happen to need very high roots, like the 200th root of 10, then this is still an interesting technique. On the other end, for roots of numbers less than 3, all the approximations are very good. Since this is not limited to integers, should you need to generate roots of, say, 1.2, this technique will give extremely close approximations.
I did realize that the convergence can be greatly improved by raising the order of x in both the numerator and denominator, effectively creating additional terms which only have an effect for small x. But, since I started this exploration with the observation that there was a sequence which was easy to remember and compute in my head, I decided to stay with that theme and stop here, rather than create much more complicated expressions. (At this point, were any of the non-math-interested folks still reading, they would go "What, you mean the previous equations weren't complicated?!" :-)
This is where I was last night, but this morning I realized that it could be a MacLaurin series. A little browsing in the CRC book and sure enough, it is the series for ey, with y=a⁄b. So, the actual limit is ea⁄b, or solving for a⁄b, we want a⁄b = ln N; indeed, ln 2 is around .693, and 9/13 is a good approximation for it. Now that I know exactly what irrational number I am trying to have for the ratio a⁄b to converge to any given N, I can use continued fractions to choose good a and b to approximate that ratio.
By the way, to get an exact result at x=1, c should be set to b - a⁄N-1.
However, using Excel to plot the actual values as a function of x, it becomes clear that the convergence becomes slow for even N around 10. (7, 10, 3) gives values within .015 of 2, but the first few terms when converging to 10 are as low as 8! One can mitigate this somewhat by raising c such that the values for x=1 and x=2 are poor (do we really need an approximation to the 1st root?), bringing the errors for the next terms down, but this only reduces the error by about half.
So, if you happen to need very high roots, like the 200th root of 10, then this is still an interesting technique. On the other end, for roots of numbers less than 3, all the approximations are very good. Since this is not limited to integers, should you need to generate roots of, say, 1.2, this technique will give extremely close approximations.
I did realize that the convergence can be greatly improved by raising the order of x in both the numerator and denominator, effectively creating additional terms which only have an effect for small x. But, since I started this exploration with the observation that there was a sequence which was easy to remember and compute in my head, I decided to stay with that theme and stop here, rather than create much more complicated expressions. (At this point, were any of the non-math-interested folks still reading, they would go "What, you mean the previous equations weren't complicated?!" :-)
Tuesday, June 30, 2009
Rational Approximation to Roots
I've been playing with a new math observation. It started when I noticed the following sequence for the nth roots of 2: (1 + 7⁄17)2 ≈ 2;
(1 + 7⁄27)3 ≈ 2; (1 + 7⁄37)4 ≈ 2; (1 + 7⁄47)5 ≈ 2. This can be generalized to (1 + 7⁄(10x-3))x = ((10x + 4)⁄(10x-3))x ≈ 2.
So, is this an isolated fluke, or are there expressions like this for all roots? It turns out there are multiple expressions (1 + a⁄(bx-c))x which converge to any N for large x. For 2, the combinations (2, 3, 1), (5, 7, 2), (7, 10, 3) and (9, 13, 4) give increasingly closer convergence to 2. The ratio a⁄b seems to be about .694, but I haven't yet mastered the limit expressions to get a symbolic expression for the number. I've found through trial and error combinations for roots of 3, 4, 5, and 10 - however, the convergence seems slower for larger N, meaning that the approximations are very poor for the lower roots (those you're more likely to use).
So, I haven't yet determined if this discovery is significant or trivial, but if you should happen to be at the sort of party where people are impressed by mental math (I haven't yet found one), you can rattle off that the 7th root of 2 is roughly 1 + 9⁄87. Or, more likely, you'll just have a handy small-number ratio for a root if you're doing math by hand.
(1 + 7⁄27)3 ≈ 2; (1 + 7⁄37)4 ≈ 2; (1 + 7⁄47)5 ≈ 2. This can be generalized to (1 + 7⁄(10x-3))x = ((10x + 4)⁄(10x-3))x ≈ 2.
So, is this an isolated fluke, or are there expressions like this for all roots? It turns out there are multiple expressions (1 + a⁄(bx-c))x which converge to any N for large x. For 2, the combinations (2, 3, 1), (5, 7, 2), (7, 10, 3) and (9, 13, 4) give increasingly closer convergence to 2. The ratio a⁄b seems to be about .694, but I haven't yet mastered the limit expressions to get a symbolic expression for the number. I've found through trial and error combinations for roots of 3, 4, 5, and 10 - however, the convergence seems slower for larger N, meaning that the approximations are very poor for the lower roots (those you're more likely to use).
So, I haven't yet determined if this discovery is significant or trivial, but if you should happen to be at the sort of party where people are impressed by mental math (I haven't yet found one), you can rattle off that the 7th root of 2 is roughly 1 + 9⁄87. Or, more likely, you'll just have a handy small-number ratio for a root if you're doing math by hand.
Subscribe to:
Posts (Atom)
Posts
