Enhanced polynomial support

[GoogleCodeExporter]

I propose to reverse the order of polynomial coefficients used by polyval
etc. It was changed to the present order for compatibility with Matlab and
SciPy, but I find this convention annoying every time I use it.

First, the Matlab absolutely makes no sense when working with polynomials
as truncations of power series. Second, having entry i be the coefficient
of x^i and nothing else is just more convenient and less error-prone.

Original issue reported on code.google.com by fredrik....@gmail.com on 18 Nov 2008 at 5:46

[GoogleCodeExporter]

I really have to say that I'm very used to this ordering:

cn*x^n + ... c1*x + c0

I'm even so used to it, that it would be more error-prone to me if it was 
inversed.

Original comment by Vinzent.Steinberg@gmail.com on 18 Nov 2008 at 9:04

[GoogleCodeExporter]

That's the ordering that *isn't* used currently. Right now it's c0*x^n + ... + 
c{n-
1}*x + cn.

Original comment by fredrik....@gmail.com on 18 Nov 2008 at 9:13

[GoogleCodeExporter]

Sorry, I meant exactly this ordering. I did not name the coefficients correctly.

Original comment by Vinzent.Steinberg@gmail.com on 18 Nov 2008 at 9:18

[GoogleCodeExporter]

I'm also used to reading polynomials that way. But it's very cumbersome for 
indexing.

Original comment by fredrik....@gmail.com on 18 Nov 2008 at 9:22

[GoogleCodeExporter]

Why? For sparse polynomials?

Original comment by Vinzent.Steinberg@gmail.com on 18 Nov 2008 at 9:47

[GoogleCodeExporter]

It's cumbersome because many algorithms involve explicitly indexing the 
coefficient
of x^k, which with the Matlab convention becomes N-k instead of just k. The 
Matlab
convention is especially annoying when building a polynomial by repeatedly 
appending
higher order terms, since each new term changes the indexing of all preceding 
terms.

But how about implementing a class for 1D polynomials? It would be fairly easy 
to
write a parser so you could do p = poly('3.25*x^4 + 5*x + 1'), which would 
obviously
be much less error prone than [1,5,0,0,3.25] or [3.25,0,0,5,1].

Original comment by fredrik....@gmail.com on 16 Mar 2009 at 6:30

[GoogleCodeExporter]

That's a nice idea and easy to implement. We could even have 
poly.append_high_order()
and poly.append_low_order() or something to solve the problems you described.

What about (re-)using sympy's polys?

Original comment by Vinzent.Steinberg@gmail.com on 16 Mar 2009 at 6:51

[GoogleCodeExporter]

Implementing a polynomial class from scratch is not much work. I've done it at 
least
five times (and interestingly, it turns out entirely different every time :-).

This is a stab at the skeleton code. The features should be self-explanatory 
(most is
documented). Some issues, apart from coefficient order, are whether polys 
should be
mutable (probably -- not much need for hashing them), whether .coeffs should 
return a
copy or give a reference to the in-place list, whether __len__ and __getitem__ 
should
be implemented, etc.

Switching to a sparse representation could perhaps be useful. An interesting 
problem
would be to generalize the Horner code to handle sparse polynomials optimally 
(there
is a faster way to evaluate x + x^100 than to perform 100 multiplications).

I added a polyfit function that essentially the same as Matlab's. However, it 
strikes
me that there is no particular need to restrict a function like this to 
polynomials
-- it could be made more general, to work with any user-supplied basis functions
(polynomials by default), and the interface could be merged with that for 
Chebyshev
approximation.

Original comment by fredrik....@gmail.com on 18 Mar 2009 at 4:37

• Changed title: Enhanced polynomial support

Attachments:

• polynomials.py

[GoogleCodeExporter]

Critical missing features: parser, adapting other functions (taylor, chebyfit, 
etc)
to work with polys, and edit affected tests/doctests.

Original comment by fredrik....@gmail.com on 18 Mar 2009 at 4:48

[GoogleCodeExporter]

Nice!

How do you want to fit nonlinear functions? Using overdetermined equation 
systems?
If it can't be formulated as a linear system, we'll need a nonlinear solver with
global convergence.

Original comment by Vinzent.Steinberg@gmail.com on 18 Mar 2009 at 7:25

[GoogleCodeExporter]

The polyfit algorithm works for any linear combination of user-supplied basis
functions, C1*f1(x)+C2*f2(x)+.... Perhaps also for a quotient of polynomials, 
or more
generally a quotient of linear combinations (C1*f1(x)+...)/(D1*g1(x)+...)
(multiplying by the denominator to obtain a linear system).

In general, nonlinear fitting indeed amounts to general nonlinear optimization.

Original comment by fredrik....@gmail.com on 18 Mar 2009 at 7:33

[GoogleCodeExporter]

Why the string magic in horner_code()? Do you think it's faster than a loop?

Original comment by Vinzent.Steinberg@gmail.com on 19 Mar 2009 at 5:01

[GoogleCodeExporter]

Yes. Also, it could be used to generate C code (at least when the coefficients 
are
real) without any changes.

Original comment by fredrik....@gmail.com on 19 Mar 2009 at 5:05

[GoogleCodeExporter]

>>> poly(range(100)).compile()
s_push: parser stack overflow
------------------------------------------------------------
Traceback (most recent call last):
  File "<ipython console>", line 1, in <module>
  File "polynomials.py", line 182, in compile
    return horner_fun(self.coeffs[::-1])
  File "polynomials.py", line 76, in horner_fun
    return eval('lambda x:'+horner_code(c))
MemoryError

However, I don't think anyone needs a polynomial of this degree.

Original comment by Vinzent.Steinberg@gmail.com on 19 Mar 2009 at 5:17

[GoogleCodeExporter]

Actually, I get a MemoryError already with range(34) which is not that big. It 
could
certainly be useful to take more than 33 terms of a Taylor series, for example.

In any case, there should be more options. For example, it should be possible to
compile a polynomial to a fixed-point version for fast high precision 
evaluation.

Original comment by fredrik....@gmail.com on 19 Mar 2009 at 5:31

[GoogleCodeExporter]

[deleted comment]

[GoogleCodeExporter]

In [26]: c = map(mpf, range(10))

In [27]: f2 = lambda x: horner(c, x)

In [28]: f1 = poly(c).__call__

In [29]: %timeit f2(7)
1000 loops, best of 3: 304 µs per loop

In [30]: %timeit f1(7)
1000 loops, best of 3: 318 µs per loop

Original comment by Vinzent.Steinberg@gmail.com on 19 Mar 2009 at 6:24

• Changed state: Started

Attachments:

• horner.py