np.linalg.cond gives incorrect results for large ill-conditioned matrices (Trac #1846)

[numpy-gitbot]

Original ticket http://projects.scipy.org/numpy/ticket/1846 on 2011-05-27 by trac user rgrout, assigned to @pv.

Hi,

I notice the following behaviour.

In [31]: np.linalg.cond(np.vander(np.arange(1,16)))
Warning: invalid value encountered in power
Warning: invalid value encountered in power
Warning: invalid value encountered in power
Warning: invalid value encountered in power
Warning: invalid value encountered in power
Warning: invalid value encountered in power
Warning: invalid value encountered in power
Out[31]: 1166991362111.4319

Where as Matlab and Octave agree on

octave:2> cond(vander(1:15))
ans =  2.5824e+21

Look at the source code, both numpy and matlab/octave attempt to find the condition number in the same way (by calculating the SVD). However, the SVD that numpy returns is completely different than matlab/octave's SVD.

When I give np.vander a small range (np.arange(10)), numpy returns a condition number that agrees with octave.

[numpy-gitbot]

@charris wrote on 2011-05-27

Works for me in the development branch, so should work in 1.6 also. This is using ATLAS and 64 bit linux.

In [5]: np.linalg.cond(np.vander(np.arange(1,16)))
Out[5]: 2.5824105551466631e+21

[numpy-gitbot]

@pv wrote on 2011-05-27

Works for me:

>>> import numpy as np
>>> np.linalg.cond(np.vander(np.arange(1,16)))
2.5824105551466631e+21

The fact that you get strange results from SVD probably implies that the linear algebra libraries (LAPACK and BLAS) you have installed are somehow faulty.

So: what is your platform, and which versions of LAPACK and BLAS do you have installed? I'm with LAPACK 3.3.0 and BLAS is provided by atlas 3.8.3, on 64-bit Linux.

[numpy-gitbot]

trac user rgrout wrote on 2011-05-27

I built numpy 1.6 from scratch and still receive incorrect results

In [6]: np.linalg.cond(np.vander(np.arange(1,16)))
/usr/lib/python2.7/site-packages/numpy/lib/twodim_base.py:518: RuntimeWarning: invalid value encountered in power
  X[:,i] = x**(N - i - 1)
Out[6]: 1166991362111.4319

My platform is Linux 2.6.38-ARCH #1 SMP PREEMPT i686 Intel(R) Pentium(R) M processor 1600MHz GenuineIntel GNU/Linux

BLAS: 3.3.0 LAPACK: 3.3.0 I do not have atlas installed

[numpy-gitbot]

@charris wrote on 2011-05-27

OK, you have 32 bit longs and the vander matrix is generated with integer powers of integers, so you have overflow generating negative integers. On my machine longs are 64 bit.

I think it might be considered an error for vander to use integers.

In [28]: vander(arange(1,16)).dtype
Out[28]: dtype('int64')

Try np.arange(1,16, dtype=float). You can also use polyvander:

In [16]: from numpy.polynomial.polynomial import polyvander

In [17]: np.linalg.cond(polyvander(np.arange(1,16), 14))
Out[17]: 5.8958225106645792e+20

Which automatically converts the arguments to floats, but the columns are in reverse order from vander. Note that the difference between the condition numbers is due to the columns being reversed, which tells you something about the accuracy of the result.

[numpy-gitbot]

trac user rgrout wrote on 2011-05-27

Replying to [comment:4 charris]:

OK, you have 32 bit longs and the vander matrix is generated with integer powers of integers, so you have overflow generating negative integers. On my machine longs are 64 bit.

I think it might be considered an error for vander to use integers.

{{{ In [28]: vander(arange(1,16)).dtype Out[28]: dtype('int64') }}}

Thanks for the help. Your suggestion worked perfectly.

It seems like the vander() function could be modified as follows.

x = asarray(x, dtype=float) #<- specify dtype when calling asarray
    if N is None:
        N=len(x)
    X = ones( (len(x),N), x.dtype)
    for i in range(N - 1):
        X[:,i] = x**(N - i - 1)
    return X

[numpy-gitbot]

@charris wrote on 2011-05-27

There is some talk of deprecating the current polynomials and replacing them with the polynomial package now in numpy. But we should probably keep vander for matlab compatibility, so a fix like that looks like the right thing. If you are looking to do large degree polynomial fits, I'd recommend the Chebyshev polynomials in the polynomial package.

Actually, I'm not exactly clear as to why the error is being raised, positive integer powers of negative integers should be OK and I don't think overflow is being checked. Errors in casting a float to integer perhaps?