Bug in Graph.is_chordal

[Jan]

Hi,

The Graph.is_chordal function can be used to return a chordless cycle
on non-chordal graphs, but the implementation seems to be wrong.

Example (ran on the sagenb.org server):
sage: g = Graph({3:[2,1,4],2:[1],4:[1],5:[2,1,4]})
sage: g.is_chordal()
=> False
sage: _, g1 = g.is_chordal(certificate=True); g1.is_chordal()
=> True

Currently, the function proceeds as follows:
1. run algorithm LexBFS to find a vertex ordering which would be a
perfect-elimination-order (PEO) if G is chordal.
2. check if this ordering satisfies the PEO-property w.r.t. G, and if
not, remember the last vertex (i) in the PEO which violates it.
(assume G is not chordal)
3. let G' be the vertex induced subgraph of G corresponding to all
vertices occuring after (i) in the PEO.
4. find two non-adj. vertices (j), (k) in G', which were both adj. to
(i) in G.
5. find a shortest path between them in G', and return the vertex
induced subgraph consisting of this path plus (i).

The resulting cycle is not necessarily chordless, since the shortest
path can contain another neighbour (x) of (i), resulting in the chord
{x,i}.
One linear-time solution is given in [1]. I think it is not hard to
modify the code to use this.

[1] Addendum: Simple Linear-Time Algorithms to Test Chordality of
Graphs, Test Acyclicity of Hypergraphs, and Selectively Reduce Acyclic
Hypergraphs.
http://dx.doi.org/10.1137/0214020

Regards,
Jan

[Nathann Cohen]

Hello Jan !

I just wrote a patch for that. What they do in this paper is exactly the same, except that of course they compute the shortest path in the graph without the neighbors of v, short of the two we are interested in :-)

http://trac.sagemath.org/sage_trac/ticket/11735

The patch is now waiting for a review. I also added a small test in the code, so that the certificate is checked for correctness before being returned. It takes absolutely no time to check, so... :-)

Oh, and thank you very much for this bug report ! Could I know what you are doing with chordal graphs ?

Nathann

[Jan]

Hi Nathann,

Thanks for the reply. As for my interest in the subject, I am an
undergrad. CS student at TU Delft and have been following a research
track there, focusing on planning/scheduling problems. We use chordal
graphs to maintain certain information related to constraints,
especially about shortest paths.

On topic: I think this does not make the procedure correct.
Also, it seems function lex_BFS does something else than what is
specified in the comments: the predecessor tree returns the _last_
discoverer for each vertex, instead of the first.
But in function is_chordal the code reads:
if t_peo.out_degree(v)>0 and g.neighbors(v) not in
neighbors_subsets[t_peo.neighbor_out_iterator(v).next()]:
it seems here we expect that the discoverer recorded is indeed the
last.

If we change the lex_BFS function to return the first discoverers,
graph({1:[2,3,4,5],2:[3,5],4:[3,5]}) is a counterexample.
I don't know whether the current impl. is correct, but i.m.o. it is
quite different from what is proposed in the paper so it might be
better to change it to something which has been proven correct.
I tried to implement the method from the paper myself, the resulting
diff is below.

Regards,
Jan

==========
9498c9498,9500
<
---
>
> pos_in_peo = dict(zip(peo, range(self.order())))
>
9500c9502
< for v in peo:
---
> for v in reversed(peo):
9502,9504c9504
< if t_peo.out_degree(v)>0 and g.neighbors(v) not in
neighbors_subsets[t_peo.neighbor_out_iterator(v).next()]:
<
< x = t_peo.neighbor_out_iterator(v).next()
---
> if t_peo.out_degree(v)>0 and [v1 for v1 in g.neighbors(v) if pos_in_peo[v1] > pos_in_peo[v]] not in neighbors_subsets[t_peo.neighbor_out_iterator(v).next()]:
9508d9507
<
9513,9519c9512,9519
<
< S = neighbors_subsets[x]
<
< for y in g.neighbors(v):
< if [y] not in S:
< break
<
---
> max_tup = (-1, 0)
> for xi in g.neighbors(v):
> for yi in g.vertices():
> if not [yi] in neighbors_subsets[xi] and pos_in_peo[xi] > pos_in_peo[v] and pos_in_peo[yi] > pos_in_peo[v]:
> new_tup = (pos_in_peo[xi], pos_in_peo[yi])
> if max_tup < new_tup:
> max_tup = new_tup
> x, y = xi, yi
9538,9539d9537
< g.delete_vertex(v)
<
==========

[Jan]

Your code, after the patch, seems to work (I tested it on all graphs
with 8 vertices, and it doesn't fail), but I think it differs from
what the paper does.
The first difference is that, after LexBFS, the current code processes
the vertices in the PEO order, and chooses the first violating vertex.
In the paper, they choose the _last_.
Also, their way of choosing the two neighbours is different.

The above code I posted is wrong (new try is below).
My new implementation passes all graphs with 8 vertices, too. This
seems hard to test ...

==========
9488c9488

<
---
> pos_in_peo = dict(zip(peo, range(self.order())))

9490,9492c9490
< for v in peo:

<
< if t_peo.out_degree(v)>0 and g.neighbors(v) not in
neighbors_subsets[t_peo.neighbor_out_iterator(v).next()]:

---
> for v in reversed(peo):

9494c9492

< x = t_peo.neighbor_out_iterator(v).next()
---
> if t_peo.out_degree(v)>0 and [v1 for v1 in g.neighbors(v) if pos_in_peo[v1] > pos_in_peo[v]] not in neighbors_subsets[t_peo.neighbor_out_iterator(v).next()]:

9503,9509c9501,9510

<
< S = neighbors_subsets[x]
<
< for y in g.neighbors(v):
< if [y] not in S:
< break
<
---
> max_tup = (-1, 0)

> nb1 = [u for u in g.neighbors(v) if pos_in_peo[u] > pos_in_peo[v]]
> for xi in nb1:
> for yi in nb1:
> if not [yi] in neighbors_subsets[xi]:

> new_tup = (pos_in_peo[xi], pos_in_peo[yi])
> if max_tup < new_tup:
> max_tup = new_tup
> x, y = xi, yi

> g.delete_vertices([vv for vv in g.vertices() if pos_in_peo[vv] < pos_in_peo[v]])
9528,9529d9528
< g.delete_vertex(v)
<
==========

[Nathann Cohen]

Hello Jan !!!

 

I am sorry for my veeeery slow answers, I am on vacation right now, with very very bad WiFi connections when I get some. If you think you would sleep better with copying the implementation given in the paper, then the best is probably to write a patch for this. I like mine better, just because I feel like I understand how it works, but to be honest I do not really mind in this case as it is sooooo easy to check whether the code is correct. What would you think of writing a patch to change the current behaviour to match the paper using your code, while letting my version of it (the updated/fixed one) in the code as a comment ?

 

Before returning the result, it should be checked, for instance like I did in my patch. If at some point the value returned is incorrect, the exception should be

"There was an error in the computed answer... Please report the bug" or something alike so that we quickly learn of it and fix the mistake. What do you think ?

 

What I fear the most are silent mistakes.. The ones you do not notice. In these situations I don't mind Sage answering me that it wasn't able to do the computations, especially when the patch is already written :-)

 

If you have time to create the patch, I will try to review it as soon as I get back to the civilisation (possibly next friday). Otherwise we'll talk about it then :-)

 

Thank you for your work, by the way !

 

Nathann

[Nathann Cohen]

Hello !

I'm finally back to the civilisation if you want us to deal with this patch :-)

Nathann

[Jan]

Hi Nathann,

after this line (line 9520 in http://trac.sagemath.org/sage_trac/attachment/ticket/11735/trac_11735.patch):
g.delete_vertices([vv for vv in g.neighbors(v) if vv != y and vv !=
x])
How can one be sure that x and y are still connected? (otherwise no
path x -- y exists)
I don't know whether this holds for the current code.
In the paper another method of choosing v,x,y is described, and the
authors prove that x,y stay connected.
For the patch, I'm not familiar with the version control system used
and unfortunately I don't have time to learn it now.

Regards,
Jan

On Sep 1, 6:23 pm, Nathann Cohen <nathann.co...@gmail.com> wrote:
> Hello !
>
> I'm finally back to the civilisation if you want us to deal with this patch :-)
>
> Nathann
>

[Nathann Cohen]

Hellooooooo Jan !

It is 20:20, it is almost dark outside and I am getting home by bike.
Worst of all, I am being assaulted by hungry mosquitoes. I re-read it
the following enough times to be sure that if I miswrote "v" instead
of "x" somewhere, reading it again wouldn't have changed anything.

When v is considered for removal, it is a leaf of the lex-BFS tree.
Its father (and first discoverer) is named x, and we suppose that
there is a vertex y which is not a neighbor of x, otherwise v is
removed without any problem.
Why should there be a x-y path avoiding all of v's neighbors ? Well,
first note that v is the furthest vertex from the source. Not
strictly, of course, but a lex-BFS is a BFS, and as we are testing for
removal the vertices in the reverse of the discovery order, it means
that in the lex-BFS exploration the vertex v is the last one to be
discovered (without considering the vertices that have been deleted
since the beginning of the elimination algorithm). Hence, it is the
one which is the further away from the source, even if other vertices
may be at equal distance from the source.
Now, how can v be adjacent to the other vertices of the graph ? Well,
for instance v can not be connected to any of x's ancestors, as its
distance from the source is strictly greater than the distance between
x and the source. Hence there is a path from x to the source of the
lex-BFS which does not touch any of x's neighbors. By the same reason,
the vertex y can not be closer than v from the source. If it were,
there would be a path from y to the source (from ancestor to ancestor)
which would avoid all of v's neighbors.

Hence, at some point, when the lex-BFS algorithm was processing all
the vertices at distance d(source, y) = d(source, v), it picked y
before v (we are removing the vertices in reverse order). But then, we
know that y and v have different sets of neighbors among the vertices
at distance "d(source, { v or y} ) - 1" from the source. We know that
because by assumption y is not adjacent to x. In this case, because y
was picked before v in the lex-BFS ordering, it means that the lex-BFS
code of y at this moment was strictly greater than the lex-BFS code of
v.

Hence there is a neighbor of y at distance d(source, y) - 1 which is
not a neighbor of v. Hence there is a path from y to the source going
through this vertex which avoids all of v's neighbors, hence
connectedness, hence certificate, hence I can jump on my bike to
escape the mosquitoes.

Nathann

[Jan]

Hi Nathann,

First, thank you for taking time to give a very detailed reply.
I'm sorry but I'm not yet done bothering you :-]
I think this is wrong:

> When v is considered for removal, it is a leaf of the lex-BFS tree.
> Its father (and first discoverer) is named x, and we suppose that
> there is a vertex y which is not a neighbor of x, otherwise v is
> removed without any problem.

(I also mentioned this in message #3, but it seems it still holds)
Consider: g = graph({1:[2,3,4,5],2:[3,5],4:[3,5]}).
g is like [the 4-cycle 2--3--4--5] + [(1,x) for x in [2 .. 5]]
g is not chordal.
A LexBFS ordering (reverse elimination ordering) could be [1, 2, 3, 5,
4].
Notice that all v in [2 .. 5] have father equal to 1.
But then x is adjacent to every other vertex, i.e. we can't find y.

Regards,
Jan

[Nathann Cohen]

> Consider: g = graph({1:[2,3,4,5],2:[3,5],4:[3,5]}).

Hahahahhahaaha ! Dead right :-)

And the code works anyway because the tree it returns actually is *NOT* a BFS tree :-)

sage: g.lex_BFS(tree = True)[1].edges()
[(2, 1, None), (3, 2, None), (4, 5, None), (5, 2, None)]

Two (obvious) black patches in my previous proof :
* It totally ignores the edges between classes of vertices with the same distance to the source -- precisely what your counter-example destroys.
* without considering "the last neighbor of v discovered before v itself", its process does not ensure that the neighborhood of a removed vertex is indeed simplicial. Very bad indeed :-)

For a time I thought I would be able to patch my proof, but in the end I don't get how this recognition algorithm works... Which is a pity because I am sure I had found a clear explanation in a paper when I implemented all that stuff. What I need right now is some sleep though. Sounds like this will really require a patch :-)

Thank you for your perseverance ! At the least I enjoyed the time thinking of this algorithm again and (re ?)-learned a few nice tricks ;-)

Good niiiiiiiiight !

Nathann

[Jan]

Hi Nathann,

> For a time I thought I would be able to patch my proof, but in the end I
> don't get how this recognition algorithm works... Which is a pity because I
> am sure I had found a clear explanation in a paper when I implemented all
> that stuff. What I need right now is some sleep though. Sounds like this
> will really require a patch :-)

I found it quite difficult to find a paper or text which explains the
chordless cycle finding.
After a long search I got only one (see post #1), but thats probably
not the one you are talking about.
So the current code with patch #11735 appears to work, but it seems
better to patch it, especially since the lexbfs function is also
patched ( #11736 http://trac.sagemath.org/sage_trac/ticket/11736 ).
As I mentioned, I don't have time to learn the patching process right
now, but I made a diff of generic_graph.py in post # 4.
I hope it's not so much work to generate a patch from that.

I would be very interested if you again find the reference which
explains how to find the cycle, it must be really hard to find (or I'm
really bad at finding references -_-).

Regards,
Jan

[Nathann Cohen]

Hello Jan !

Quick update : 

I just reviewed #11738 which touched some code related to chordality.

http://trac.sagemath.org/sage_trac/ticket/11738

I also updated my patch at #11735 that implements my (still unproved) version of the algorithm, because it is formally better than the current one, and most importantly because it checks that the results returned are *CORRECT*, which is urgent indeed. At least if this algorithm does not work either, an exception will be raised instead of returning a wrong result.

I will produce a patch matching your code/paper as soon as possible, though as I also have to prepare my defense's talk and move to belgium before the beginning of November, I cannot do it just now `:-D`

I do swear this will be fixed however, and that *before* I begin to work in Brussels. But I will already sleep easier if I know Sage doesn't return wrong results ! Thank you again for your help ! I'll keep this thread updated about what's happening with this code.

Nathann

[Nathann Cohen]

Hello Jan !

I am trying to write the patch corresponding to your modifications,
but the file has changed much since and I have some trouble dealing
with your diff file... Could you copy/paste the totality of the code ?
:-)

Thaaaaaaaaaanks !

Nathann

[Jan]

Hi Nathann,

I extracted the modified is_chordal function, it is shown below.

Regards, Jan

==========
def is_chordal(self, certificate = False):
r"""
Tests whether the given graph is chordal.

A Graph `G` is said to be chordal if it contains no induced hole
(a
cycle of length at least 4).

Alternatively, chordality can be defined using a Perfect
Elimination
Order :

A Perfect Elimination Order of a graph `G` is an ordering
`v_1,...,v_n`
of its vertex set such that for all `i`, the neighbors of `v_i`
whose
index is greater that `i` induce a complete subgraph in `G`.
Hence, the
graph `G` can be totally erased by successively removing vertices
whose
neighborhood is a clique (also called *simplicial* vertices)
[Fulkerson65]_.

(It can be seen that if `G` contains an induced hole, then it can
not
have a perfect elimination order. Indeed, if we write
`h_1,...,h_k` the
`k` vertices of such a hole, then the first of those vertices to
be
removed would have two non-adjacent neighbors in the graph.)

A Graph is then chordal if and only if it has a Perfect
Elimination
Order.

INPUT:

- ``certificate`` (boolean) -- Whether to return a certificate.

* If ``certificate = False`` (default), returns ``True`` or
``False`` accordingly.

* If ``certificate = True``, returns :

* ``(True, peo)`` when the graph is chordal, where ``peo``
is a
perfect elimination order of its vertices.

* ``(False, Hole)`` when the graph is not chordal, where
``Hole`` (a ``Graph`` object) is an induced subgraph of
``self`` isomorphic to a hole.

ALGORITHM:

This algorithm works through computing a Lex BFS on the graph,
then
checking whether the order is a Perfect Elimination Order by
computing
for each vertex `v` the subgraph induces by its non-deleted
neighbors,
then testing whether this graph is complete.

This problem can be solved in `O(m)` [Rose75]_ ( where `m` is the
number
of edges in the graph ) but this implementation is not linear
because of
the complexity of Lex BFS. Improving Lex BFS to linear complexity
would
make this algorithm linear.

The complexity of this algorithm is equal to the complexity of the
implementation of Lex BFS.

EXAMPLES:

The lexicographic product of a Path and a Complete Graph
is chordal ::

sage: g =
graphs.PathGraph(5).lexicographic_product(graphs.CompleteGraph(3))
sage: g.is_chordal()
True

The same goes with the product of a random lobster
( which is a tree ) and a Complete Graph ::

sage: g =
graphs.RandomLobster(10,.5,.5).lexicographic_product(graphs.CompleteGraph(3))
sage: g.is_chordal()
True

The disjoint union of chordal graphs is still chordal::

sage: (2*g).is_chordal()
True

Let us check the certificate given by Sage is indeed a perfect
elimintion order::

sage: (_, peo) = g.is_chordal(certificate = True)
sage: for v in peo:
... if not g.subgraph(g.neighbors(v)).is_clique():
... print "This should never happen !"
... g.delete_vertex(v)
sage: print "Everything is fine !"
Everything is fine !

Of course, the Petersen Graph is not chordal as it has girth 5 ::

sage: g = graphs.PetersenGraph()
sage: g.girth()
5
sage: g.is_chordal()
False

We can even obtain such a cycle as a certificate ::

sage: (_, hole) = g.is_chordal(certificate = True)
sage: hole
Subgraph of (Petersen graph): Graph on 5 vertices
sage: hole.is_isomorphic(graphs.CycleGraph(5))
True

REFERENCES:

.. [Rose75] Rose, D.J. and Tarjan, R.E.,
Algorithmic aspects of vertex elimination,
Proceedings of seventh annual ACM symposium on Theory of
computing
Page 254, ACM 1975

.. [Fulkerson65] Fulkerson, D.R. and Gross, OA
Incidence matrices and interval graphs
Pacific J. Math 1965
Vol. 15, number 3, pages 835--855

TESTS:

Trac Ticket #11735::

sage: g = Graph({3:[2,1,4],2:[1],4:[1],5:[2,1,4]})

sage: _, g1 = g.is_chordal(certificate=True); g1.is_chordal()

False
sage: g1.is_isomorphic(graphs.CycleGraph(g1.order()))
True
"""

# If the graph is not connected, we are computing the result on
each component
if not self.is_connected():

# If the user wants a certificate, we had no choice but to
# collect the perfect elimination orders... But we return
# a hole immediately if we find any !
if certificate:
peo = []
for gg in self.connected_components_subgraphs():

b, certif = gg.is_chordal(certificate = True)
if not b:
return certif
else:
peo.extend(certif)

return True, peo

# One line if no certificate is requested
else:
return all( gg.is_chordal() for gg in
self.connected_components_subgraphs() )

peo,t_peo = self.lex_BFS(reverse=True, tree=True)

g = self.copy()

# Remembering the (closed) neighborhoods of each vertex
from sage.combinat.subset import Subsets
neighbors_subsets = dict([(v,Subsets(self.neighbors(v)+[v])) for v
in self.vertex_iterator()])
pos_in_peo = dict(zip(peo, range(self.order())))
# Iteratively removing vertices and checking everything is fine.
for v in reversed(peo):

if t_peo.out_degree(v)>0 and [v1 for v1 in g.neighbors(v) if
pos_in_peo[v1] > pos_in_peo[v]] not in
neighbors_subsets[t_peo.neighbor_out_iterator(v).next()]:

# Do we need to return a hole ?
if certificate:

# In this case, let us take two nonadjacent neighbors
of
# v. In order to do so, we pick a vertex y which is a
# neighbor of v but is not adjacent to x, which we
know
# exists by the test written two lines above.

max_tup = (-1, 0)
nb1 = [u for u in g.neighbors(v) if pos_in_peo[u] >
pos_in_peo[v]]
for xi in nb1:
for yi in nb1:
if not [yi] in neighbors_subsets[xi]:
new_tup = (pos_in_peo[xi], pos_in_peo[yi])
if max_tup < new_tup:
max_tup = new_tup
x, y = xi, yi
g.delete_vertices([vv for vv in g.vertices() if
pos_in_peo[vv] < pos_in_peo[v]])

g.delete_vertices([vv for vv in g.neighbors(v) if vv !
= y and vv != x])

g.delete_vertex(v)

# Our hole is v + (a shortest path between x and y not
# containing v or any of its neighbors).

hole = self.subgraph([v] + g.shortest_path(x,y))

# There was a bug there once, so it's better to check
the
# answer is valid, especally when it is so cheap ;-)

if hole.order() <= 3 or not hole.is_regular(k=2):
raise Exception("The graph is not chordal, and
something went wrong in the computation of the certificate. Please
report this bug, providing the graph if possible !")

return (False, hole)
else:
return False

if certificate:
return True, peo

else:
return True

[Nathann Cohen]

Hellooooooooooo again !!!

I finally wrote this patch, available there :

http://trac.sagemath.org/sage_trac/ticket/11961

Though the two algorithms are very similar, I thought it better to
split it into two different parts, so that one can understand better
how each of them works without having to bear with "if algorithm == A"
every second line. Two implementations are now available, two are
checked before being returned, both pass all tests. We should be on
the safe side, now :-D

Thank you for your help !! :-)

Nathann

[Jan]

Hi Nathann,

Thanks for fixing it! :) I forgot the issue so my reply is late.

I looked at the patch you produced.
I think it's fine except for that I would prefer the "algorithm"
keyword-argument to be placed at the end of the argument list for
backward compatibility.

Best regards,
Jan

[Nathann Cohen]

Hello !!

> Thanks for fixing it! :) I forgot the issue so my reply is late.

Well, it took me a month to write it ;-)

> I looked at the patch you produced.
> I think it's fine except for that I would prefer the "algorithm"
> keyword-argument to be placed at the end of the argument list for
> backward compatibility.

Right ! I just updated the patch. Thank you again !

Nathann