Minimax Polynomial Fit



#3

Hi all,

Busy as always, just a quick post to let you know that a new article
from mine appears published in the latest September-October issue of Datafile (V24 N5), namely:

HP-71B Minimax Polynomial Fit

This is a 12-page article including two programs, namely a 3-liner to compute the
Least-Squares Nth-degree polynomial which best fits a given set of (x,y) datapoints
in the least-square sense, and a 50-line program to compute the minimax Nth-degree
polynomial fit to an arbitrary set of datapoints (x,y), which can be entered from the
keyboard, automatically generated by evaluating an arbitrary user-provided function,
or read/stored from/to a data file.

The user can specify the degree of the polynomial or else specify instead a maximum
absolute error and let the program select the lowest degree which will achieve it.
Once the coefficients have been computed, the user can evaluate the resulting
polynomial at any given arguments right from the keyboard (even as part of more
complicated expressions), and can also scan any given interval for maximum errors in
it, at whatever step size desired.

For example, when asked to automatically produce the minimax polynomial of lowest
degree which fits y = Gamma(x) in the interval [1,2] with maximum absolute error less
than 0.00005 (i.e.: 4 correct decimal places), my program
produces the 5th-degree minimax polynomial:

    P(x) = A0 + A1*x + A2*x^2 + ... + A5*x^5
with
    A0 =  3.73459956
A1 = -6.72738793
A2 = 6.54439953
A3 = -3.37095678
A4 = 0.91923039
A5 = -0.09992158
Once computed, we can use the automatic scan option to check the [1,2] interval at 0.01 steps
to confirm that it has absolute maximum error app. 0.00003812 (<0.00005) at x=1.06
    >FNE(1, 2, 0.01)

1.06000000 (x-coordinate where the max. err. occurs)
0.00003812 (absolute maximum error in [1,2])

and can evaluate it for any argument, say Ln(Pi), like this:

    >FIX 5 @ FNF(LN(PI))

0.93479 (P(Ln(Pi)) = predicted value of Gamma(Ln(Pi)))

while the correct value is
    GAMMA(LN(PI)) -> 0.93480
i.e., we sure got 4 correct decimal places, more nearly five.

Any and all comments are very welcome.

Best regards from V.


#4

Quote:
The user can specify the degree of the polynomial or else specify instead a maximum absolute error and let the program select the lowest degree which will achieve it.

Hello Valentin,

Thanks for letting us know of your new article. Just one question:
can the user also specify polynomials with odd or even powers only? This would be particularly useful when fitting odd and even functions, as less coefficients would be needed for a given accuracy.

Best regards,

Gerson.


#5

Hi, Gerson:

Gerson posted:

"Can the user also specify polynomials with odd or even powers only? This would be particularly
useful when fitting odd and even functions, as less coefficients would be needed for a given accuracy."

First of all, thanks for your interest in my article. As for your question, there's no need of a separate option to achieve that effect. You just simply specify an interval where it shows that the function is odd or even, and the minimax error condition itself will guarantee that only odd or even powers are used.

For instance, let's say that you want to fit a minimax polynomial of lowest possible degree to y=sin(x) for values of x not exceeding Pi/2 and using just odd powers of x (as sin(-x) = -sin(x), so sin(x) is an odd function). Also, you want a guaranteed maximum error <= 1E-6.

As I've stated, you only need to specify an interval where the 'oddness' of the functions shows, in this case [-Pi/2, Pi/2] will
do. For this particular function, interval, and maximum error,
my program returns the following 7th-degree polynomial:

P(x) = 0.99999661 * x - 0.16664831 * x3 + 0.00830636 * x5
- 0.00018365 * x7

with guaranteed absolute maximum error less than 1E-6 in [-Pi/2, Pi/2]. As you can see, only odd powers of x were returned, all the remaining coefficients being zero within machine rounding limits so our minimax approximation only requires four non-zero coefficients.

Here's a short table of sin(x) vs. P(x):

>FIX 7 @ FOR X=-PI/2 TO PI/2 STEP PI/20 @ DISP X,SIN(X),FNF(X),RES-SIN(X) @ NEXT X

X SIN(X) P(X) ERROR
-------------------------------------------------------------------------------
-1.5707963 -1.0000000 -0.9999994 0.0000006
-1.4137167 -0.9876883 -0.9876887 -0.0000004
-1.2566371 -0.9510565 -0.9510560 0.0000005
-1.0995574 -0.8910065 -0.8910061 0.0000004
-0.9424778 -0.8090170 -0.8090173 -0.0000003
-0.7853982 -0.7071068 -0.7071074 -0.0000006
-0.6283185 -0.5877853 -0.5877856 -0.0000003
-0.4712389 -0.4539905 -0.4539903 0.0000002
-0.3141593 -0.3090170 -0.3090164 0.0000006
-0.1570796 -0.1564345 -0.1564340 0.0000005
0.0000000 0.0000000 0.0000000 0.0000000
0.1570796 0.1564345 0.1564340 -0.0000005
0.3141593 0.3090170 0.3090164 -0.0000006
0.4712389 0.4539905 0.4539903 -0.0000002
0.6283185 0.5877853 0.5877856 0.0000003
0.7853982 0.7071068 0.7071074 0.0000006
0.9424778 0.8090170 0.8090173 0.0000003
1.0995574 0.8910065 0.8910061 -0.0000004
1.2566371 0.9510565 0.9510560 -0.0000005
1.4137167 0.9876883 0.9876887 0.0000004
1.5707963 1.0000000 0.9999994 -0.0000006

and, as you can see, the absolute error never exceeds 0.0000006 < 1E-6, so the results are accurate to six decimal places plus or minus 1 ulp.

This exact same approach, but specifying the [-Pi/6, Pi/6] interval instead, yields a 7th-degree, odd-powers-only, minimax polynomial which has a guaranteed absolute maximum error e = 3.2E-11 in that interval, while using only four coefficients.

Combined with range reduction, this could form the basis of a simple yet very accurate implementation of trigonometric functions in machines lacking them, as you know. Similar approaches would allow implementing most any specialized function, say Gamma(x), Lambert's W(x), etc.

Best regards from V.

Edited: 11 Oct 2005, 5:55 a.m.


#6

Quote:
As I've stated, you only need to specify an interval where the 'oddness' of the functions shows, in this case [-Pi/2, Pi/2] will do. For this particular function, interval, and maximum error, my program returns the following 7th-degree polynomial:


P(x) = 0.99999661 * x - 0.16664831 * x3 + 0.00830636 * x5 - 0.00018365 * x7


Hi Valentin,

Thanks for your detailed and comprehensive reply. However, I still have another question, if you don't mind: is it possible to force the first coefficient to 1 and then recalculate the other coefficients? Just rounding the first coefficient to 1 would cause the accuracy to decrease, as we can see in the table below:

      P(x) = x - 0.16664831 * x3 + 0.00830636 * x5 - 0.00018365 * x7


    X                   SIN(X)               P(X)                  ERROR
-------------------------------------------------------------------------------
0.0000000 0.0000000 0.0000000 0.0000000
0.1570796 0.1564345 0.1564345 0.0000001
0.3141593 0.3090170 0.3090175 0.0000005
0.4712389 0.4539905 0.4539919 0.0000014
0.6283185 0.5877853 0.5877877 0.0000024
0.7853982 0.7071068 0.7071100 0.0000032
0.9424778 0.8090170 0.8090205 0.0000035
1.0995574 0.8910065 0.8910099 0.0000033
1.2566371 0.9510565 0.9510602 0.0000037
1.4137167 0.9876883 0.9876935 0.0000051
1.5707963 1.0000000 1.0000046 0.0000046

Of course the new polynomial with recalculated coefficients wouldn't be as accurate as the original one, but would be better than this.
Setting the first coefficient to 1 is important to reduce the number of steps or registers necessary to implement these polynomials on a RPN calculator.

Best regards,

Gerson.


#7

Hi again, Gerson:

Gerson wrote:
"is it possible to force the first coefficient to 1 and then recalculate the other coefficients?"

Yes, there are two ways of doing it:

  • Modifying the program itself: as one of the coefficients would be fixed, you'll actually need one equation less to solve, and so one less datapoint in your tuples. A very simple ad-hoc modification to the SPOLMM subprogram, slightly affecting the part where it is building the matrix of the system to solve, and then including the fixed coefficient in the returned coefficient's array would suffice. This would be faster and most accurate, but requires modifying the program as it is now.

  • You can read the reference you give in your own article, which I did, where it is explictily stated how to proceed to achieve exactly what you want, even giving a fully worked example to make it crystal-clear.

To save you the trouble of reading it, basically you transform the problem of finding a 7th-degree minimax polynomial for y=sin(x) in [0, Pi/4], say, for the related problem of finding a 2nd-degree polynomial for y=f(x) in [0, (Pi/4)^2], where f(x) is defined thus:

    f(x) = sin(sqrt(x))/x^(3/2)-1/x
this way you get a polynomial which directly translates to your desired special-form one for the original y=sin(x) function in the original [0, Pi/4] interval.

You can achieve that with my program by simply selecting [D]efined function, then answering the "f(x)=" prompt like this:

    f(x) = SIN(SQR(X))/(X*SQR(X))-1/X
and specifying the [0, (Pi/4)^2] interval, as stated. You'll get the results you want.

Note:As f(x) is defined, it would generate an error when evaluated at x=0. You can either specify the interval using a small number (1E-6, say) instead of 0, or else include a user-defined function, like this:

    DEF FNY(X) @ IF X=0 THEN FNY=-1/6 ELSE FNY=SIN(SQR(X))/(X*SQR(X))-1/X
then answer the "f(x)=" prompt with f(x)=FNY(X)

Best regards from V.

Edited: 13 Oct 2005, 6:37 a.m.


#8

Hi Valentin,

Thanks for your step by step explanations, they've really helped me. And thanks for reminding me to read my own reference too :-)
At the time I didn't give too much importance to these details because I didn't have any minimax tool available. As I said, I used curve fitting and did some minor adjustments by hand to get the best results I could. I got a 9.58E-11 maximum error with the three coefficients I obtained for SIN this way. The coefficients I obtained using a true minimax polynomial yielded a maximum 9.76E-11 error. Not bad! (Of course, my manual solution took much more time and would be impossible to be done for higher degree polynomials).

Using the procedure you showed me, I found the odd power coefficients for a 13th degree polynomial to approximate ATAN (of course, the coefficient indeces refer to the auxiliary 5th degree polynomial):

8 DEF FNY(Y) @ IF X=0 THEN FNY=-1/3 ELSE FNY=ATAN(SQR(X))/(X*SQR(X))-1/X
9 END DEF

X1,X2,N=0,3-2*SQR(2),20
.
.
.
H= 0.00000000050
A0=-0.33333333283
A1= 0.19999978075
A2=-0.14284193118
A3= 0.11072111863
A4=-0.08629967887
A5= 0.05061006653

Just for a comparison, here are the same coefficients calculated by an expensive mathematical software, using the tecnique shown in the reference:
a0=-0.3333333328040560
a1= 0.1999997758846710
a2=-0.1428416640124080
a3= 0.1107161126831530
a4=-0.0862630685349892
a5= 0.0505192312223211

The slight difference observed is due, of course, to different machine precisions. Anyway, the maximum absolute errors are 5E-10 and 7E-11, respectively. A remarkable achievement for both programmer and HP-71B! You program is a joy to use, I wish I only had a real 71B!

Best regards,

Gerson.

P.S.: Below are my approximations to ATAN and SIN. Maximum errors
observed in a spread sheet were 6.6E-11 and 9.8E-11 respectively,
but I am curious to see how they would behave in a HP-12CP. Any
simulator available?

---------------------------------------------------------------------


a = -0.3333333329
b = 0.1999997763
c = -0.142841662
d = 0.1107161
e = -0.0862632
f = 0.05052

atan(x) = x*(1+x^2*(a+x^2*(b+x^2*(c+x^2*(d+x^2*(e+f*x^2)))))),

[0 <= x <= sqrt(2)-1]

atan(a) = 2*atan(sqr(1/(a^2)+1)-1/a)

---------------------------------------------------------------------

a = -0.1666666652
b = 0.008333217
c = -0.000197282


sin(x) = x*(1+x^2*(a+x^2*(b+c*x^2))), [0 <= x <= pi/6]

sin(a) = sin(a/3)*(3-4*sin^2(a/3))


---------------------------------------------------------------------


#9

Hi, Gerson:

Gerson posted:
"The slight difference observed is due, of course, to different machine precisions. Anyway, the maximum absolute errors are 5E-10 and 7E-11, respectively. A
remarkable achievement for both programmer and HP-71B! You program is a joy to use, I wish I only had a real 71B!"


I'm glad you liked my program and found it useful. Owning a physical HP-71 is a real pleasure but certainly this kind of program benefits a lot from being run in an emulator (say Jean-François Garnier's Emu71 or Hrastprogrammer's HP-71X), because it's so math intensive. This way you get results in seconds and can try different approaches with utmost ease, making the 'fitting experience' rather enjoyable, IMHO.

Also, it's a real testimony to HP-71B/Math capabilities that this kind of functionality, which usually requires advanced software to accomplish, can be implemented in a mere 50 lines of code, and moreover, providing three different input methods (including mass
storage) and two modes of output. A basic implementation would surely fit in 30 lines of code or so.

And it's minimaxing, no less ! :-)

Best regards from V.

#10

Valentin --

Looks like an impressive effort! However, I'll have to delve into my HP-71B manuals in order to try out your work.

I saw no mention of this program needing a Math ROM, but I have one with the manual. Now all I need is a HP-71B Surveying Pac manual to go with the ROM...

-- KS


#11

Hi, Karl:

Karl posted:

"I saw no mention of this program needing a Math ROM, but I have one with the manual."

Thanks for your continued interest in my humble productions. Yes, I didn't mention that a Math ROM is needed, but of course it is.

If you've followed my posts since long, you'll know that when I write "HP-71B" it goes without saying that it includes a Math ROM, because if not, it is *not* the real HP-71B originally designed by its wise creators, but the maimed thingy they were forced to release when it was clear the System ROMs were going to exceed the allotted 64 Kb of ROM space.

The original HP-71B did include most functionality now present in the Math ROM, including complex numbers definition and functionality, of which some scraps still remain in a bare-bones HP-71B (try to enter the program line: 10 A=(2,3)+SIN((4,5)) and see it being duly accepted, no syntax complaints), and so for me an HP-71B deserves that name if and only if it includes a Math ROM.

Else, you should change its name, perhaps "HP-71BL" for "Light" will do, though "HP-71BM" would be better, or even "HP-71BFU" ... ;-)

BTW, don't miss my upcoming S&SCM #11, due today ! :-)

Best regards from V.


Possibly Related Threads...
Thread Author Replies Views Last Post
  Polynomial program Michael Carey 12 2,316 11-20-2013, 08:43 PM
Last Post: CR Haeger
  Best statistical fit Richard Berler 8 1,692 10-30-2013, 11:25 PM
Last Post: Walter B
  HP Prime polynomial long division bluesun08 13 2,336 10-30-2013, 03:29 AM
Last Post: parisse
  [HP-Prime xCAS] Review Polynomial Tools + BUGs + Request CompSystems 0 605 09-05-2013, 12:53 PM
Last Post: CompSystems
  Matrix Characteristic Polynomial - Reloaded. Ángel Martin 12 2,121 08-22-2013, 05:33 PM
Last Post: Thomas Klemm
  WP34s program submission: Quadratic fit Andrew Nikitin 2 812 06-13-2013, 02:44 AM
Last Post: Paul Dale
  New Blog Post: Length of a Polynomial Segment Eddie W. Shore 1 765 01-17-2013, 08:56 PM
Last Post: Namir
  New empirical fit for ln(x) Namir 0 507 12-11-2012, 12:49 PM
Last Post: Namir
  Where does the 17B-II fit in? Matt Agajanian 22 3,564 03-25-2012, 10:52 AM
Last Post: Don Shepherd
  Where does the 38G fit in? Matt Agajanian 5 1,078 03-10-2012, 01:08 PM
Last Post: Marcus von Cube, Germany

Forum Jump: