[SciPy-user] Characteristic polynomial

Tue Dec 11 11:38:17 EST 2007

Hello Bryan and thanks for responding.
I have already tried that, but check this example...

>>> from numpy import *
>>> a=arange(25)
>>> a.shape=5,5
>>> a
array([[ 0,  1,  2,  3,  4],
       [ 5,  6,  7,  8,  9],
       [10, 11, 12, 13, 14],
       [15, 16, 17, 18, 19],
       [20, 21, 22, 23, 24]])
>>> coeffs=poly(a)
>>> coeffs
array([  1.00000000e+00,  -6.00000000e+01,  -2.50000000e+02,
        -1.67542241e-13,   1.82969746e-28,  -8.57128973e-45])
>>> poly1d(coeffs)
poly1d([  1.00000000e+00,  -6.00000000e+01,  -2.50000000e+02,
        -1.67542241e-13,   1.82969746e-28,  -8.57128973e-45])

Check these coefficients, like -1.67542241e-13 etc... these numbers are
supposed to be zero, but they are "very small".

The very same example in sage :
A = matrix(5,5, range(25)); A
A.charpoly()

x^5 - 60*x^4 - 250*x^3

I tried to apply round function on the coefficients :
>>> for i in range(len(coeffs)): coeffs[i]=round(coeffs[i])
...
>>> coeffs
array([   1.,  -60., -250.,    0.,    0.,    0.])

and while in this case it returns the right result, there are cases where
the coefficients are like : 1.1234566789e64,
where the results cannot be rounded (since they are already integers) but
they are different where they were supposed to be the same.
A strange example :

For the 10x10 matrices :
A=matrix([[36,8,8,1,1,1,1,1,8,1],
             [8,36,8,8,1,1,1,1,1,1],
             [8,8,25,8,1,8,1,1,1,1],
             [1,8,8,36,8,1,1,1,1,1],
             [1,1,1,8,36,1,8,8,1,1],
             [1,1,8,1,1,36,8,1,1,8],
             [1,1,1,1,8,8,49,1,1,1],
             [1,1,1,1,8,1,1,49,8,1],
             [8,1,1,1,1,1,1,8,36,8],
             [1,1,1,1,1,8,1,1,8,49]])

B=matrix([[36,8,8,1,1,1,1,8,1,1],
             [8,25,8,1,8,8,1,1,1,1],
             [8,8,36,8,1,1,1,1,1,1],
             [1,1,8,36,8,1,1,1,1,8],
             [1,8,1,8,36,8,1,1,1,1],
             [1,8,1,1,8,36,8,1,1,1],
             [1,1,1,1,1,8,49,1,8,1],
             [8,1,1,1,1,1,1,49,8,1],
             [1,1,1,1,1,1,8,8,36,8],
             [1,1,1,8,1,1,1,1,8,49]])

>>> poly(A)
array([  1.00000000e+00,  -3.88000000e+02,   6.65430000e+04,
        -6.63946000e+06,   4.26558031e+08,  -1.84256082e+10,
         5.41543001e+11,  -1.06843860e+13,   1.35292436e+14,
        -9.91717494e+14,   3.19104663e+15])
>>> poly(B)
array([  1.00000000e+00,  -3.88000000e+02,   6.65430000e+04,
        -6.63946000e+06,   4.26558031e+08,  -1.84256082e+10,
         5.41543001e+11,  -1.06843860e+13,   1.35292436e+14,
        -9.91717494e+14,   3.19104663e+15])

sage returns :

x^10 - 388*x^9 + 66543*x^8 - 6639460*x^7 + 426558031*x^6 -
18425608218*x^5 + 541543001174*x^4 - 10684385967388*x^3 +
135292435663763*x^2 - 991717494497240*x + 3191046625350150

for both matrices.

You can see the difference between 3.19104663e+15  and 3191046625350150
On the other hand, if we get the exact value of 3.19104663e+15 (by using :
poly(A)[-1] and poly(B)[-1]) we have :
3191046625350152.0
3191046625350153.5

And with all these... :S I'm really confused.
Anyway, it would be really useful if there existed a function to compute the
modular characteristic polynomial (modulo a big prime) (like in maple).

On Dec 11, 2007 5:37 PM, Bryan Van de Ven <bryanv at enthought.com> wrote:

> In [16]: numpy.lib.poly?
> Type:           function
> Base Class:     <type 'function'>
> String Form:    <function poly at 0xb7a66224>
> Namespace:      Interactive
> File:
> /workspace/python2.5/lib/python2.5/site-packages/numpy/lib/polynomial.py
> Definition:     numpy.lib.poly(seq_of_zeros)
> Docstring:
>     Return a sequence representing a polynomial given a sequence of roots.
>
>     If the input is a matrix, return the characteristic polynomial.
>
>     Example:
>
>         >>> b = roots([1,3,1,5,6])
>         >>> poly(b)
>         array([ 1.,  3.,  1.,  5.,  6.])
>
>
> In [17]: a = array([[1,0],[0,2]])
>
> In [18]: coeffs = numpy.lib.poly(a)
>
> In [19]: numpy.lib.poly1d(coeffs)
> Out[19]: poly1d([ 1., -3.,  2.])
>
>
>
>
> Nikolas Karalis wrote:
> > Hello... I'm trying to compute the characteristic polynomial of an
> > integer (numpy) matrix. But i cannot find any way of doing this.
> > Do you have any ideas/suggestions?
> >
> > --
> > Nikolas Karalis
> > Applied Mathematics and Physics Undergraduate
> > National Technical University of Athens, Greece
> > http://users.ntua.gr/ge04042
> >
> >
> > ------------------------------------------------------------------------
> >
> > _______________________________________________
> > SciPy-user mailing list
> > SciPy-user at scipy.org
> > http://projects.scipy.org/mailman/listinfo/scipy-user
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>

-- 
Nikolas Karalis
Applied Mathematics and Physics Undergraduate
National Technical University of Athens, Greece
http://users.ntua.gr/ge04042
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20071211/6055dbfb/attachment.html>