[SciPy-user] Characteristic polynomial

[Anne Archibald]

On 23/12/2007, Nikolas Karalis <nikolaskaralis at gmail.com> wrote:

> >>> A=matrix([[290, 324, 323, 364, 340, 341, 365, 336, 342, 326],
>         [324, 290, 322, 338, 366, 341, 365, 336, 344, 324],
>         [323, 322, 286, 337, 337, 366, 364, 361, 320, 323],
>         [345, 327, 326, 315, 343, 344, 342, 363, 343, 312],
>         [329, 347, 325, 343, 319, 345, 343, 364, 327, 327],
>         [329, 329, 347, 344, 345, 323, 344, 342, 348, 328],
>         [347, 347, 346, 342, 343, 344, 321, 341, 326, 311],
>         [325, 325, 343, 363, 364, 342, 341, 313, 324, 309],
>         [342, 344, 320, 362, 338, 366, 338, 335, 286, 307],
>         [340, 339, 338, 330, 361, 362, 333, 329, 316, 265]])
>
> >>> B=matrix([[290, 324, 336, 365, 341, 340, 364, 323, 342, 326],
>         [324, 290, 336, 339, 367, 340, 364, 322, 344, 324],
>         [325, 325, 313, 341, 342, 364, 363, 343, 324, 309],
>         [346, 328, 341, 319, 346, 343, 342, 347, 345, 312],
>         [330, 348, 342, 346, 321, 345, 344, 346, 329, 327],
>         [328, 328, 364, 343, 345, 321, 341, 326, 346, 328],
>         [346, 346, 363, 342, 344, 341, 317, 325, 324, 311],
>         [323, 322, 361, 365, 365, 338, 336, 286, 320, 323],
>         [342, 344, 335, 364, 340, 364, 336, 320, 286, 307],
>         [340, 339, 329, 332, 363, 360, 331, 338, 316, 265]])
>
> >>> poly(A)
> array([  1.00000000e+00,  -3.00800000e+03,  -1.10612600e+06,
>         -1.55679754e+08,  -1.10343405e+10,  -4.15092661e+11,
>          -7.89507268e+12,  -6.51631023e+13,  -1.80794492e+14,
>         -2.21111889e+00,  -2.95926167e-13])
>
> >>> poly(B)
> array([  1.00000000e+00,  -3.00800000e+03,  -1.10612600e+06,
>         -1.55679754e+08 ,  -1.10343405e+10,  -4.15092661e+11,
>         -7.89507268e+12,  -6.51631023e+13,  -1.80794492e+14,
>          4.59224482e+00,   3.39308585e-14])
>
> As you can see, the 2 polynomials are the same (up to the last 2 terms). The
> last term is considered to be "almost" 0, so we can see that the difference
> is the coefficient of x^1.
>
> If we compute the exact same thing with Maple and Sage, we get (for both
> matrices) :
> x^10 - 3008*x^9 - 1106126*x^8 - 155679754*x^7 - 11034340454*x^6 -
> 415092661064*x^5 - 7895072675601*x^4 - 65163102265268*x^3 -
> 180794492489124*x^2
>
> so, it is the same, since it doesn't "compute" the x^1 term.
>
> This also happens for a few other matrices...
>
> Can anybody help me with this? What is the right answer for the char. poly
> (i guess it is the Sage's and Maple's one, since they agree).
> Why does this difference occur? Is it "sage" to ignore the last 2 terms?

When you do all but the simplest matrix computations with numpy, it
treats the matrices as  containing floating-point numbers. Since
floating-point numbers are as a rule only approximations to the exact
value you would like ot be working with, the answers you get out are
only approximations. If you look at the coefficients you're getting
from numpy's characteristic polynomial calculation, you will see that
the coefficient of x^1 is about 10^-14 times the coefficient of x^2. A
good rule of thumb is that you can consider as almost zero any number
that is many orders of magnitude less than the biggest number in your
answer. So 4, or 2, or whatever, is "close to zero" in this context.

More specifically, when you compute a small number as the difference
of two very large numbers, that number has far fewer digits of
accuracy than the original two numbers. Many matrix computations, and
in particular eigenvalue and characteristic polynomial calculations,
involve just this kind of delicate cancellation, often many times
over. So you have to expect errors to accumulate.

Some problems, like finding the roots of a polynomial whose
coefficients you know, are just inherently "unstable", that is, if you
change the coefficients by one in the last significant figure, the
roots move around wildly. These problems are basically hopeless, and
you will need to reformulate your problem in a way that doesn't
require you to caluclate the roots. Matrix inversion can be like this,
for example, so rather than solve Ax=b as x=A^-1 b, you should use a
specialized solver (of which there are many). For your problem, the
answer depends on what you are doing with the characteristic
polynomials. If you are (for example) trying to find all the
eigenvalues, use numpy's eigenvalue solvers directly. If you only want
certain coefficients - the x^0 and x^(n-1) are the easiest - think
about whether there's a better formula (e.g. determinant or trace) for
them.

Your case is a bit special: your input numbers are not floating-point.
They are integers, and (I assume) they represent the input values
*exactly*. Thus in principle it should be possible to compute the
characteristic polynomial exactly. Unfortunately, using floating-point
is almost certainly a bad way to do this. (MAPLE might allow you to
use floating-point with (say) several hundred decimal places, which
might result in very nearly integral coefficients, but this is going
to be *slow*.) Instead, you should use a special-purpose integer
linear algebra routine. These integer linear algebra routines are
specialized, computationally expensive, and not part of numpy. They
are part of SAGE and MAPLE. More generally, this kind of exact
computation belongs to the domain of "symbolic computation", which is
not something numpy is for.

In summary: use SAGE, or think about your problem more.

Anne