[Numpy-discussion] Problem in LinearAlgebra?

[Nadav Horesh]

I tested the problem with:
     1. Numeric 23.1 under python 2.3.2
     2. numarray 0.8 (I made a copy of the Scientific package where all
        calls to Numeric were replaced to numarray), under python 2.3.2
There results where about the same -- high coefficients for the 5th
order polynomials.
I would expect reliable fit for a high order polynomials only under very
special circumstances, so this is not a big surprise. My advice is:
      * Make sure that this is a bug and not a result of a numerical
        instability. If you can trace it down and point to a bug, then
        report it. The numarray package is very usable and is under a
        very active and rapid development, thus bugs are being fixed
        fast.
      * Look for a solution in the scipy package: It is generally better
        then Scientific.
      * Polynomials fit is relatively very simple --- you may write one
        of you own in less then a one day work. Since, as I said, the
        problem is, in many cases, unstable, you'll have the chance to
        implement more stable linear-equation solvers.

Nadav.


On Fri, 2003-10-31 at 15:19, Rob W.W. Hooft wrote:
> I am using Polynomial.py from Scientific Python 2.1, together with 
> Numeric 17.1.2. This has always served me well, but now we are busy 
> upgrading our software, and I am currently porting some code to 
> Scientific Python 2.4.1, Numeric 22.0. Suddenly I do no longer manage to 
> get proper 2D polynomial fits over 4x4th order. At 5x5 the coefficients 
> that come back from LinearAlgebra.linear_least_squares have exploded. In 
> the old setup, I easily managed 9x9th order if I needed to, but most of 
> the time I'd stop at 6x6th order. Would anyone have any idea how this 
> difference can come about? I managed to work around this for the moment 
> by using the equivalent code in the fitPolynomial routine that uses 
> LinearAlgebra.generalized_inverse (and it doesn't even have problems 
> with the same data at 8x8), but this definitely feels not right! I can't 
> remember reading anything like this here before.
> 
> Together with Konrad Hinsen, I came to the conclusion that the problem 
> is not in Scientific Python, so it must be the underlying LinearAlgebra 
> code that changed between releases 17 and 22.
> 
> I hacked up a simplified example. Not sure whether it is the most simple 
> case, but this resembles what I have in my code, and I'm quite sure it 
> worked with Numeric 17.x, but currently it is horrible over order (4,4):
> 
> --------------------------------------
> import Numeric
> 
> def func(x,y):
>      return x+0.1*x**2+0.01*x**4+0.002*x**6+0.03*x*y+0.001*x**4*y**5
> 
> x=[]
> y=[]
> z=[]
> for dx in Numeric.arange(0,1,0.01):
>      for dy in Numeric.arange(0,1,0.01):
>          x.append(dx)
>          y.append(dy)
>          z.append(func(dx,dy))
> 
> from Scientific.Functions import Polynomial
> data=Numeric.transpose([x,y])
> z=Numeric.array(z)
> for i in range(10):
>      print data[i],z[i]
> pol=Polynomial.fitPolynomial((4,4),data,z)
> print pol.coeff
> ------------------------------------
> for 4,4 this prints:
> [[  1.84845529e-05  -7.60502772e-13   2.71314749e-12  -3.66731796e-12
>           1.66977148e-12]
>   [  9.99422967e-01   3.00000000e-02  -3.26346097e-11   4.42406519e-11
>          -2.01549767e-11]
>   [  1.03899464e-01  -3.19668064e-11   1.14721790e-10  -1.55489826e-10
>           7.08425891e-11]
>   [ -9.40275000e-03   4.28456838e-11  -1.53705205e-10   2.08279772e-10
>          -9.48840470e-11]
>   [  1.80352695e-02  -1.10999843e-04   8.00662570e-04  -2.17266676e-03
>           2.47500004e-03]]
> 
> for 5,5:
> 
> [[ -2.25705839e+03   6.69051337e+02  -6.60470163e+03   6.66572425e+03
>          -8.67897022e+02   1.83974866e+03]
>   [ -2.58646837e+02  -2.46554689e+03   1.15965805e+03   7.01089888e+03
>          -2.11395436e+03   2.10884815e+03]
>   [  3.93307499e+03   4.34484805e+02  -4.84080392e+03   5.90375330e+03
>           1.16798049e+03  -4.14163933e+03]
>   [  1.62814750e+03   2.08717457e+03   1.15870693e+03  -3.37838057e+03
>           3.49821689e+03   5.80572585e+03]
>   [  4.54127557e+02  -1.56645524e+03   4.58997025e+00   1.69772635e+03
>          -1.37751039e+03  -7.59726558e+02]
>   [  2.37878239e+03   9.43032094e+02   8.58518644e+02  -8.35846339e+03
>          -5.55845668e+02   1.87502761e+03]]
> 
> Which is clearly wrong.
> 
> I appreciate any help!
> 
> Regards,
> 
> Rob