On Sat, Feb 27, 2010 at 10:02 PM, Friedrich Romstedt <friedrichromst...@gmail.com> wrote: > 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'm sorry this comment turns out to be confusing. It has apparently quite the contrary effect of what I wanted to achieve: Since there is already a polynomial module in numpy I wanted to highlight their difference and stress that they are used to do arithmetic, e.g. compute f([x],[y]) = [x] * (sin([x])**2 + [y]) in Taylor arithmetic. > >> 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? no, the vectorization is something different. It's purpose becomes only clear when applied in Algorithmic Differentiation. E.g. if you have a function f: R^N -> R x -> y=f(x) where x = [x1,...,xN] and you want to compute the gradient g(x) of f(x), then you can compute df(x)/dxn by propagating the following array of Taylor polynomials: x = numpy.array( UTPS([x1_0, 0]), ..., UTPS([xn_0, 1]), ..., UTPS([xN_0,0]), dtype=object) y = f(x) if you want to have the complete gradient, you will have to repeat N times. Each time for the same zero'th coefficients [x1,...,xN]. Using the vecorized version, you would do only one propagation x = numpy.array( UTPS([x1_0, 1,0,...,0]), ..., UTPS([xn_0, 0,...,1,...0]), ..., UTPS([xN_0,0,....,1]), dtype=object) y = f(x) i.e. you save the overhead of calling the same function N times. > >> [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. Yes, I'm sorry, this is a typo. > > 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). the goal is to have several Taylor polynomials evaluated in the same base point, e.g. [x_0, x_{11}, x_{21}, x_{31}] [x_0, x_{12}, x_{22}, x_{32}] [x_0, x_{13}, x_{23}, x_{33}] i.e. P=3, D=3 One could use an (P,D) array. However, one would do unnecessary computations since x_0 is the same for all P polynomials. I.e. one implements the data structure as [x]_{D,P} := [x_0, x_{1,1},...,x_{1,D-1},x_{2,1},...,x_{P,D-1}]. This results in a non-const stride access. > > 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? I believe the degree D is typically much to small (i.e. D <= 4) to justify the additional overhead of using FFT, though there may be cases when really high order polynomials are used. > > Friedrich > _______________________________________________ > NumPy-Discussion mailing list > NumPy-Discussion@scipy.org > http://mail.scipy.org/mailman/listinfo/numpy-discussion > _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion