[Numpy-discussion] Problem in LinearAlgebra?

[David M. Cooke]

On Fri, Oct 31, 2003 at 02:19:54PM +0100, 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.

Works for me:

(4,4)
[[ 0.00001848 -0.          0.         -0.          0.        ]
 [ 0.99942297  0.03       -0.          0.         -0.        ]
 [ 0.10389946 -0.          0.         -0.          0.        ]
 [-0.00940275  0.         -0.          0.         -0.        ]
 [ 0.01803527 -0.000111    0.00080066 -0.00217267  0.002475  ]]

(5,5)
[[-0.00000175 -0.          0.          0.         -0.          0.        ]
 [ 1.00008231  0.03       -0.          0.         -0.          0.        ]
 [ 0.09914353 -0.          0.         -0.          0.         -0.        ]
 [ 0.00350289  0.         -0.          0.         -0.          0.        ]
 [ 0.00333036 -0.          0.         -0.          0.          0.001     ]
 [ 0.00594     0.         -0.          0.         -0.          0.        ]]

(6,6)
[[ 0.    -0.     0.    -0.     0.    -0.     0.   ]
 [ 1.     0.03  -0.     0.    -0.     0.    -0.   ]
 [ 0.1   -0.     0.    -0.     0.    -0.     0.   ]
 [-0.     0.    -0.     0.    -0.     0.    -0.   ]
 [ 0.01  -0.     0.    -0.     0.     0.001  0.   ]
 [-0.     0.    -0.     0.    -0.     0.    -0.   ]
 [ 0.002 -0.     0.    -0.     0.    -0.     0.   ]]

(I've set sys.float_output_suppress_small to True to get a better picture.)

I'm using Numeric 23.1, ScientificPython 2.4.3, and Python 2.3.2, from Debian
unstable. Numeric is compiled to use the system LAPACK libraries (using ATLAS
for the BLAS). This is on both Athlon and PowerPC chips.

How did you compile Numeric? With or without the system LAPACK libraries?

-- 
|>|\/|<
/--------------------------------------------------------------------------\
|David M. Cooke                      http://arbutus.physics.mcmaster.ca/dmc/
|cookedm at physics.mcmaster.ca