Pari stack growning out of bounds

[Jernej Azarija]

Consider the following program that computes the spectrum and chromatic number of a graph:

for g in graphs.nauty_geng(str(7)):
    s = g.spectrum()
    g.chromatic_number()

This works quickly and like a charm. Now consider the following program that computes something related to the spectrum and chromatic number:

for g in graphs.nauty_geng(str(7)):
    s = g.spectrum()
    sp = sum([el**2 for el in s if el > 0])
    sm = sum([el**2 for el in s if el < 0])
   
    if sm != 0 and not (sp/sm+1 <= g.chromatic_number()):
        print g.graph6_string()

The program gets stuck in computing something for a long time and in the end it dies trying to allocate a ~7GB pari stack:

  File "minw.py", line 9, in <module>
    if sm != _sage_const_0  and not (sp/sm+_sage_const_1  <= g.chromatic_number()):
  File "element.pyx", line 902, in sage.structure.element.Element.__richcmp__ (sage/structure/element.c:8480)
  File "element.pyx", line 847, in sage.structure.element.Element._richcmp (sage/structure/element.c:7930)
  File "element.pyx", line 829, in sage.structure.element.Element._richcmp_ (sage/structure/element.c:7659)
  File "element.pyx", line 874, in sage.structure.element.Element._richcmp (sage/structure/element.c:8342)
  File "/home/azi/sage-5.3.rc1/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 3755, in __cmp__
    rcmp = cmp(self.real(), other.real())
  File "/home/azi/sage-5.3.rc1/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 4395, in __cmp__
    return self._sub_(other).sign()
  File "/home/azi/sage-5.3.rc1/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 4611, in sign
    return self.sign()
  File "/home/azi/sage-5.3.rc1/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 4614, in sign
    self.exactify()
  File "/home/azi/sage-5.3.rc1/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 3466, in exactify
    self._set_descr(self._descr.exactify())
  File "/home/azi/sage-5.3.rc1/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 7591, in exactify
    left.exactify()
  File "/home/azi/sage-5.3.rc1/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 3466, in exactify
    self._set_descr(self._descr.exactify())
  File "/home/azi/sage-5.3.rc1/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 7334, in exactify
    arg.exactify()
  File "/home/azi/sage-5.3.rc1/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 3466, in exactify
    self._set_descr(self._descr.exactify())
  File "/home/azi/sage-5.3.rc1/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 7591, in exactify
    left.exactify()
  File "/home/azi/sage-5.3.rc1/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 3466, in exactify
    self._set_descr(self._descr.exactify())
  File "/home/azi/sage-5.3.rc1/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 7591, in exactify
    left.exactify()
  File "/home/azi/sage-5.3.rc1/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 3466, in exactify
    self._set_descr(self._descr.exactify())
  File "/home/azi/sage-5.3.rc1/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 7593, in exactify
    gen = left._exact_field().union(right._exact_field())
  File "/home/azi/sage-5.3.rc1/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 2276, in union
    newpol, self_pol, k = pari_nf.rnfequation(my_factor, 1)
  File "gen.pyx", line 10412, in sage.libs.pari.gen._pari_trap (sage/libs/pari/gen.c:54794)
  File "gen.pyx", line 9718, in sage.libs.pari.gen.PariInstance.allocatemem (sage/libs/pari/gen.c:50859)
  File "gen.pyx", line 10233, in sage.libs.pari.gen.init_stack (sage/libs/pari/gen.c:53888)
MemoryError: Unable to allocate 65536000000 bytes memory for PARI.


What exactly is going on? From what I've checked into the sources I do not see what is causing this issue? Am I doing something stupid in the second sage program or is this some kind of Sage-Pari bug?

Best,

Jernej

[John Cremona]

A partial answer is that somehow the elements of the spectrum are
being created as elements of QQbar, and must have high degree, so pari
is (behind the scenes) trying to compare algebraic numbers of high
degree, which is hard as it involves creating number fields of high
degree. It seems computationally hard to test if an element of QQbar
is positive, which is what you are doing.

What is the degree of the elements of the spectrum which cause the problem?

Does anyone know why QQbar is being used here?

John

> --
> You received this message because you are subscribed to the Google Groups
> "sage-devel" group.
> To post to this group, send email to sage-...@googlegroups.com.
> To unsubscribe from this group, send email to
> sage-devel+...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sage-devel?hl=en.
>
>

[Rob Beezer]

Some of your eigenvales (many?) are in QQbar, the field of algebraic numbers.  Testing equality in this field is extremely expensive for some reason.

Also, sum() is sufficiently generic, that I think it is not very efficient either.

A recent patch implements eigenvalues of symmetric matrices with floating point entries.  Which will be real for adjacency matrices of your graph.  This might be "good enough" for you, and was the orginal motivation for the patch.  Can't recall if this patch has been reviewed/merged, or not.

Rob

[Rob Beezer]

On Monday, September 24, 2012 1:29:24 PM UTC-7, John Cremona wrote:

Does anyone know why QQbar is being used here?

 
For "simple" graphs, the adjaceny matrix will be symmetric, so the eigenvalues could go into AA instead of QQbar, but I would guess the problem (comparison with zero) might persist?

Where would you suggest they go?  Right now, I think these are just eigenvalues of a matrix over QQ with those outside of QQ going into QQbar (iirc).

Rob

 

[John Cremona]

It depends on what they are going to be used for. I don't work with
graphs so cannot guess! Why not have a parameter to the spectrum()
function to allow the user to choose between some floating point
approximate values, versus exact algebraic ones. And in the latter
case I don't know why using QQbar (or AA where possible) is done
instead of using a fixed number field with an embedding. I guess it
is not deliberate but because of this:

sage: M=Matrix(QQ,2,2,[0,1,2,0])
sage: M.eigenvalues()
[-1.414213562373095?, 1.414213562373095?]
sage: lam = M.eigenvalues()
sage: type(lam[0])
<class 'sage.rings.qqbar.AlgebraicNumber'>

John

> Rob

[Jeroen Demeyer]

On 2012-09-24 22:29, John Cremona wrote:
> It seems computationally hard to test if an element of QQbar
> is positive

What does it even mean for an element of QQbar to be positive?

[Dima Pasechnik]

Indeed,  in general this makes little sense without specifying an embedding. 

However, there is more to this. Indeed, say, you have f(x)=det(xI-A), and A=A^T has rational entries.

E.g. that's the case when A is the adjacency matrix of an (undirected) graph.

Then all the roots of a polynomial in Q[x] are real. And then the question like "are all roots of f positive?", "count positive roots", etc

make perfect sense.

Dima

[Dima Pasechnik]

As suggested below, it might work if you replace the line

s = g.spectrum()

with 

s = [AA(xx) for xx in g.spectrum()]

making all the eigenvalues real (as they should be, as you generate undirected graphs)

[Dima Pasechnik]

On Tuesday, 25 September 2012 15:31:09 UTC+8, Dima Pasechnik wrote:

As suggested below, it might work if you replace the line

s = g.spectrum()

with 

s = [AA(xx) for xx in g.spectrum()]

This still doesn't really work, but 

 s = [RR(xx) for xx in g.spectrum()]

does. You'd also need to weaken your inequality a bit to account for numerical noise (RR by default has 53-bit mantissa).

Doing sp/sm+1.01 instead of sp/sm+1 produces a list of 94 graphs.

How to make this robust is another question.

HTH,

Dmitrii

[Rob Beezer]

To summarize some of this discussion more carefully:

Eigenvalues of a symmetric rational matrix should already land in the algebraic field (AA), so there is no need to convert them.  Error messages appear to come from QQbar, since much of the functionality of AA is contained in the same module as QQbar.

sage: G = graphs.GridGraph([2,3])
sage: spec = G.spectrum(); spec
[2.414213562373095?, 1, 0.4142135623730951?, -0.4142135623730951?, -1, -2.414213562373095?]
sage: spec[0].parent()
Algebraic Field

Eigenvalues of rational matrices come from roots of the characteristic polynomial - for graphs this is the adjacency matrix.  If you want more control, you can get the adjacency matrix (am) and then the factored characteristic polynomial (fcp) and deal with that.

sage: G.am().fcp()
(x - 1) * (x + 1) * (x^2 - 2*x - 1) * (x^2 + 2*x - 1)

If numerical results are "good enough" you can convert the adjacency matrix to floating-point entries and get numerical eigenvalues.

sage: A = G.am().change_ring(RDF)
sage: A.eigenvalues()
[-2.41421356237, 2.41421356237, -1.0, -0.414213562373, 1.0, 0.414213562373]

But for large graphs, the above will sometimes produce (seemingly) complex roots with near-zero imaginary parts.  The 'symmetric' (or 'hermitian') keyword will use more specialized SciPy routines that yield real eigenvalues, and also find them faster.

sage: A.eigenvalues(algorithm="symmetric")
[-2.41421356237, -1.0, -0.414213562373, 0.414213562373, 1.0, 2.41421356237]

Since computing eigenvalues is equivalent to finding roots of polynomials, you can pick your poison here.

(Maybe some documentation in the graph-spectrum method about the (new) symmetric numerical algorithms is warranted?)

Rob

[Dima Pasechnik]

On Wednesday, 26 September 2012 00:54:44 UTC+8, Rob Beezer wrote:

To summarize some of this discussion more carefully:

Eigenvalues of a symmetric rational matrix should already land in the algebraic field (AA), so there is no need to convert them.  Error messages appear to come from QQbar, since much of the functionality of AA is contained in the same module as QQbar.

sage: G = graphs.GridGraph([2,3])
sage: spec = G.spectrum(); spec
[2.414213562373095?, 1, 0.4142135623730951?, -0.4142135623730951?, -1, -2.414213562373095?]
sage: spec[0].parent()
Algebraic Field

they should, but they don't. Look:

sage: a=AA(1)

Algebraic Real Field

sage: b=QQbar(1)

sage: b.parent()

Algebraic Field

 So this is a bug...

[Rob Beezer]

Thanks, Dima.

On Tuesday, September 25, 2012 10:08:13 AM UTC-7, Dima Pasechnik wrote:

 So this is a bug...

Maybe more like: my mistake. I had my print versions of the complex/real algebraic numbers mixed up. 

Here is the essence of how these (exact) eigenvalues are computed for rational matrices (not just adjacency matrices of graphs).  It would be a definite improvement to recognize symmetric matrices (through a keyword or an optional check) and then create these roots in AA rather than QQbar.

sage: G = graphs.GridGraph([3,2])
sage: p = G.am().fcp(); p

(x - 1) * (x + 1) * (x^2 - 2*x - 1) * (x^2 + 2*x - 1)

sage: f = p[2][0]; f
x^2 - 2*x - 1
sage: evs = f.root_field('t').gen(0).galois_conjugates(QQbar); evs
[-0.4142135623730951?, 2.414213562373095?]
sage: evs[0].parent()
Algebraic Field
sage: evs = f.root_field('t').gen(0).galois_conjugates(AA)
sage: evs[0].parent()
Algebraic Real Field