[Numpy-discussion] [ANN]: Taylorpoly, an implementation of vectorized Taylor polynomial operations and request for opinions

Friedrich Romstedt friedrichromstedt@gmail....
Sat Feb 27 15:02:03 CST 2010


To the core developers (of numpy.polynomial e.g.): Skip the mess and
read the last paragraph.

The other things I will post back to the list, where they belong to.
I just didn't want to have off-topic discussion there.

> I wanted to stress that one can do arithmetic on Taylor polynomials in
> a very similar was as with complex numbers.

I do not understand completely.  What is the analogy with complex
numbers which one cannot draw to real numbers?  Or, more precise: What
/is/ actually the analogy besides that there are operations?  With i,
i * i = -1, thus one has not to discard terms, contrary to the
polynomial product as you defined it, no?

> I guess there are also situations when you have polynomials   z =
> \sum_{d=0}^{D-1} z_d t^d, where z_d are complex numbers.

>> I like it more to implement operators as overloads of the __mul__
>
> I thought about something like
> def __mul__(self, other):
>    return mul(self,other)

Yes, I know!  And in fact, it may symmetrise the thing a bit.

>> etc., but this is a matter of taste.
>>  In fact, you /have/ to provide
>> external binary operators, because I guess you also want to have
>> numpy.ndarrays as left operand.  In that case, the undarray will have
>> higher precedence, and will treat your data structure as a scalar,
>> applying it to all the undarray's elements.
>
> well, actually it should treat it as a scalar since the Taylor
> polynomial is something like a real or complex number.

Maybe I misunderstood you completely, but didn't you want to implement
arrays of polynomials using numpy?  So I guess you want to apply a
vector from numpy pairwise to the polynomials in the P-object?

> [z]_{D+1} &=_{D+1}& [x]_{D+1} y_{D+1} \\
> \sum_{d=0}^D z_d t^d & =_{D+1} & \left( \sum_{d=0}^D x_d t^d \right)
> \left( \sum_{c=0}^D y_c t^c \right)

Did your forget the [] around the y in line 1, or is this intentional?
 Actually, I had to compile the things before I could imagine what you
could mean.

Why don't you use multidimensional arrays?  Has it reasons in the C
accessibility?  Now, as I see it, you implement your strides manually.
 With a multidimensional array, you could even create arrays of shape
(10, 12) of D-polynomials by storing an ndarray of shape (10, 12, D)
or (D, 10, 12).

Just because of curiosity:  Why do you set X_{0, 1} = ... X_{0, P} ?

Furthermore, I believe there is some elegant way formulating the
product for a single polynomial.  Let me think a bit ...

For one entry of [z]_E, you want to sum up all pairs:
    x_{0} y{E} + ... + x{D} y{E - D} ,     (1)
right?  In a matrix, containing all mutual products x_{i} y{j}, this
are diagonals.  Rotating the matrix by 45° counterclockwise, they are
sums in columns.  Hmmm... how to rotate a matrix by 45°?

Another fresh look: (1) looks pretty much like the discrete
convolution of x_{i} and y_{i} at argument E.  Going into Fourier
space, this becomes a product.  x_i and y_i have finite support,
because one sets x_{i} and y_{i} = 0 outside 0 <= i <= D.  The support
of (1) in the running index is at most [0, D].  The support in E is at
most [0, 2 D].  Thus you don't make mistakes by using DFT, when you
pad x and y by D zeros on the right.  My strategy is: Push the padded
versions of x and y into Fourier space, calculate their product,
transform back, cut away the D last entries, and you're done!

I guess you mainly want to parallelise the calculation instead of
improving the singleton calculation itself?  So then, even non-FFT
would incorporate 2 D explicit python sums for Fourier transform, and
again 2 D sums for transforming back, is this an improvement over the
method you are using currently?  And you could code it in pure Python
:-)

I will investigate this tonight.  I'm curious.  Irrespective of that I
have also other fish to fry ... ;-)

iirc, in fact DFT is used for efficient polynomial multiplication.
Maybe Chuck or another core developer of numpy can tell whether
numpy.polynomial does it actually that way?

Friedrich


More information about the NumPy-Discussion mailing list