infinite loop in spectral_radius for trees

[dcoudert]

Issue raised on ​https://ask.sagemath.org/question/43496/why-cant-i-find-the-spectral-radius-of-a-tree/

The following instruction will never end

sage: g = graphs.CompleteBipartiteGraph(1,3)
sage: g.spectral_radius()

Adding a print in the while loop to get e_max, e_min, c_prec, e_max-e_min, e_max * c_prec, we get

sage: g.spectral_radius()
3.0 1.0 1e-10 2.0 3e-10
3.0 1.0 1e-10 2.0 3e-10
3.0 1.0 1e-10 2.0 3e-10
...

and its the same for every tree. It's working well if the graph has a single edge, but for the path graph of order 3, we already get the infinite loop.

[dcoudert]

I don't know how to fix that.

[vdelecroix]

Indeed, this is very bad for any bipartite graph because the two eigenvalues are lambda1 (the spectral radius) and -lambda1. Something special needs to be done in that case. One simple option is to consider the square of the adjacency matrix (reduced to one of the bi partition). I will try to provide a fix.

[vdelecroix]

Though, in the case of trees as requested on ask I am sure there are much more clever things to do.

[vdelecroix]

More precisely, here is the problem. Iterating the incidence matrix gets you into a cycle of length two

sage: g = graphs.CompleteBipartiteGraph(1,3)
sage: A = g.adjacency_matrix().change_ring(RDF)
sage: v = vector([1,2,3,4])
sage: v = A*v; v/= v.norm(); v
(0.9819805060619657, 0.1091089451179962, 0.1091089451179962, 0.1091089451179962)
sage: v = A*v; v/= v.norm(); v
(0.18898223650461365, 0.566946709513841, 0.566946709513841, 0.566946709513841)
sage: v = A*v; v/= v.norm(); v
(0.9819805060619656, 0.10910894511799618, 0.10910894511799618, 0.10910894511799618)
sage: v = A*v; v/= v.norm(); v
(0.18898223650461363, 0.5669467095138409, 0.5669467095138409, 0.5669467095138409)

But with A^2 everything is fine

sage: w = A**2 * v
sage: [w[i] / v[i] for i in range(4)]
[3.0, 3.0, 3.0, 3.0]

(and the eigenvalue is indeed sqrt(3) for the initial graph).

[vdelecroix]

@dcoudert: do you know why is_bipartite is not there on directed graph? The notion of bipartiteness does not care about orientation at all.

[vdelecroix]

New commits:

​cb000cd

26148: fix spectral radius for bipartite graphs

[dcoudert]

Replying to vdelecroix:

@dcoudert: do you know why is_bipartite is not there on directed graph? The notion of bipartiteness does not care about orientation at all.

It's certainly because noone asked for this feature yet?

We can move the method to generic_graph.py. The part that tests bipartiteness is the same. However, we must update the method for finding a cycle when the digraph is not bipartite as shortest path methods care about orientation.

[dcoudert]

The test if not G treats the empty graph only, so it's not appropriate for graph without edges. The correct test is if not G.size().

The current behavior is that a graph without edges is considered bipartite. I agree with that. For the empty graph, it's always ambiguous in definitions...

Concerning the method you propose for bipartite graphs in fact you use a method similar to project_left in BipartiteGraph, except that the projection don't accept loops/multiedges.

[vdelecroix]

Replying to dcoudert:

The test if not G treats the empty graph only, so it's not appropriate for graph without edges. The correct test is if not G.size().

The current behavior is that a graph without edges is considered bipartite. I agree with that. For the empty graph, it's always ambiguous in definitions...

Concerning the method you propose for bipartite graphs in fact you use a method similar to project_left in BipartiteGraph, except that the projection don't accept loops/multiedges.

Indeed. Note that project_left also produces a graph and not a digraph. Why on earth is this called project_left?

[dcoudert]

Well, there is also method project_right ;)

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​5544dda	26148: move is_bipartite to generic_graph.py
​c4aeaca	26148: use is_bipartite when g is oriented

[vdelecroix]

I removed the call to shortest path.

[dcoudert]

• Why not writing directly return sqrt(e_min), sqrt(e_max) ?

• while s is not None: -> while s:
• u_to_root -> v_to_root (you start from v)
• before the loop while u_to_root and w_to_root and u_to_root[-1] == w_to_root[-1] you could add a comment like # Search for the nearest common ancestor of v and w. No obligation.

• You must add special case for the graph with a single vertex and no edge (don't know the expected result).

  sage: Graph(1).spectral_radius()
  ---------------------------------------------------------------------------
  ValueError                                Traceback (most recent call last)
  <ipython-input-8-c8d7e35649a7> in <module>()
  ----> 1 Graph(Integer(1)).spectral_radius()
  
  /home/dcoudert/sage/local/lib/python2.7/site-packages/sage/graphs/base/static_sparse_graph.pyx in sage.graphs.base.static_sparse_graph.spectral_radius (build/cythonized/sage/graphs/base/static_sparse_graph.cpp:12977)()
      984                         H.add_edge(u0, u2)
      985 
  --> 986         e_min, e_max = spectral_radius(H, prec)
      987         e_min = sqrt(e_min)
      988         e_max = sqrt(e_max)
  
  /home/dcoudert/sage/local/lib/python2.7/site-packages/sage/graphs/base/static_sparse_graph.pyx in sage.graphs.base.static_sparse_graph.spectral_radius (build/cythonized/sage/graphs/base/static_sparse_graph.cpp:12235)()
      957     """
      958     if not G:
  --> 959         raise ValueError("empty graph")
      960     if G.is_directed():
      961         if not G.is_strongly_connected():
  
  ValueError: empty graph
  

  Note that Graph([(0, 0)], loops=True).spectral_radius() returns (1.0, 1.0)

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​8150104	26148: renaming and comments in is_bipartite
​7c8240d	26148: fix one vertex case for spectral_radius

[vdelecroix]

Replying to dcoudert:

I did implement everything excepted

• while s is not None: -> while s:

Nope, s = 0 or s = False may be valid vertices.

[dcoudert]

Nope, s = 0 or s = False may be valid vertices.

Right :P

[dcoudert]

LGTM.