Polynomial implementation in numpy.ma
Hi all, There were some failures in the polynomial tests earlier today, and while investigating I saw that numpy.ma implements its own root finder. It uses inversion of a Van der Monde matrix, which I believe may suffer from some numerical instability problems. Given that Charles has gone to some length to implement good polynomial root finders, I think it would be best to employ those instead, and simply pre-filter the data that comes from the masked array module, if possible. Regards Stéfan
2011/9/14 Stéfan van der Walt <stefan@sun.ac.za>
Hi all,
There were some failures in the polynomial tests earlier today, and while investigating I saw that numpy.ma implements its own root finder. It uses inversion of a Van der Monde matrix, which I believe may suffer from some numerical instability problems. Given that Charles has gone to some length to implement good polynomial root finders, I think it would be best to employ those instead, and simply pre-filter the data that comes from the masked array module, if possible.
The test failure arises because the test compares against np.polyfit where Travis changed the column scaling, which in turn changed the singular values, and it is that comparison that fails. The Vandermonde matrix needs to be used for the fitting so nothing should be changed there. To fix the test failure, either the test could be modified or the ma.polyfit routine can be changed. Probably the latter is the best option, the new scaling code can be brought over and maybe the weight and covariance options added.Travis' changes needs to be audited in any case to make sure all the corner cases are covered. Chuck
On Sat, Sep 17, 2011 at 10:40 AM, Charles R Harris < charlesr.harris@gmail.com> wrote:
2011/9/14 Stéfan van der Walt <stefan@sun.ac.za>
Hi all,
There were some failures in the polynomial tests earlier today, and while investigating I saw that numpy.ma implements its own root finder. It uses inversion of a Van der Monde matrix, which I believe may suffer from some numerical instability problems. Given that Charles has gone to some length to implement good polynomial root finders, I think it would be best to employ those instead, and simply pre-filter the data that comes from the masked array module, if possible.
The test failure arises because the test compares against np.polyfit where Travis changed the column scaling, which in turn changed the singular values, and it is that comparison that fails.
The Vandermonde matrix needs to be used for the fitting so nothing should be changed there. To fix the test failure, either the test could be modified or the ma.polyfit routine can be changed. Probably the latter is the best option, the new scaling code can be brought over and maybe the weight and covariance options added.Travis' changes needs to be audited in any case to make sure all the corner cases are covered.
Note that the enhanced np.polyfit can be called directly with the weight set to zero for the masked rows, i.e., set w to the compliment of the mask. That should simplify the code a bit. Chuck
On Sat, Sep 17, 2011 at 9:40 AM, Charles R Harris <charlesr.harris@gmail.com> wrote:
The Vandermonde matrix needs to be used for the fitting so nothing should be changed there.
I remember now where I heard this: http://people.maths.ox.ac.uk/trefethen/mythstalk.pdf Specifically, the statement: "Some well-known algorithms take O(n^2) operations (e.g. naïve Lagrange formula) or are exponentially unstable (e.g. Vandermonde matrices based on monomials x^k )." I wasn't sure whether to interpret that as the Vandermonde-based algorithm being flawed, or the polynomial fitting problem on an equidistant grid being ill-conditioned in general. Regards Stéfan
2011/9/17 Stéfan van der Walt <stefan@sun.ac.za>
On Sat, Sep 17, 2011 at 9:40 AM, Charles R Harris <charlesr.harris@gmail.com> wrote:
The Vandermonde matrix needs to be used for the fitting so nothing should be changed there.
I remember now where I heard this:
http://people.maths.ox.ac.uk/trefethen/mythstalk.pdf
Specifically, the statement: "Some well-known algorithms take O(n^2) operations (e.g. naïve Lagrange formula) or are exponentially unstable (e.g. Vandermonde matrices based on monomials x^k )."
I wasn't sure whether to interpret that as the Vandermonde-based algorithm being flawed, or the polynomial fitting problem on an equidistant grid being ill-conditioned in general.
Equispaced grid points are bad for power series fits -- the Chebyshev points are better -- but power series also tend to have problems even when Chebyshev points are used since they never become particularly orthogonal over the sample points. Alternatives like Bernstein polynomials, Legendre polynomials, and Chebyshev polynomials are much better in that regard but not generally as easy to hand as power series, which one reason I wanted to make available in numpy. However, doing a fit using Chebyshev polynomials and then converting it to a power series throws away the Chebyshev advantage as the Chebyshev series is much less sensitive to numerical errors in the coefficients. Chuck
participants (2)
-
Charles R Harris
-
Stéfan van der Walt