[SciPy-user] Characteristic polynomial
Matthieu Brucher
matthieu.brucher at gmail.com
Tue Dec 11 11:45:03 EST 2007
You can't expect scipy to find the exact coefficients. sage is a symbolic
package, it will find the correct answers, but scipy will find only an
approximated one, up to the machine precision. This is what you see in your
exemple.
If you have integers, you could expect scipy to return long integers (exact
result), but this is not the case as almost everything is converted into a
float array before the actual C or Fortran routine is run.
Matthieu
2007/12/11, Nikolas Karalis <nikolaskaralis at gmail.com>:
>
> 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
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>
>
--
French PhD student
Website : http://matthieu-brucher.developpez.com/
Blogs : http://matt.eifelle.com and http://blog.developpez.com/?blog=92
LinkedIn : http://www.linkedin.com/in/matthieubrucher
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20071211/4bb02224/attachment.html>
More information about the SciPy-User
mailing list