total unimodularity testing - algorithm?

[Dima Pasechnik]

There is a total unimodularity test for matrices implemented in the matroids
mega-patch.
Where can one find a description of the algorithm used there?

Thanks,
Dima

[Rudi Pendavingh]

Hi Dima,

I don't think we discuss the implementation in much detail in the documentation strings.

In case you wonder, the algorithm we implemented is *not* the poly-time algorithm that follows from regular matroid decomposition. It is a pretty generic algorithm that also serves to test if a matroid is representable over other partial fields.

To test if the matrix that represents a matroid M is TU, we consider each flat F of rank r(M)-2. Contracting such a flat in the represented matroid M, we get a represented minor N=M/F of rank 2. Say this N is represented by the matrix
[ 1 0 1 .. 1 ..]
[ 0 1 1 .. a_i..]
If one of the a_i is not in {0,1}, then the representation of M is not TU. Conversely, if M passes this test for each F, then its representation is TU.

So it is an inherently exponential algorithm, but fairly optimized in its use of data structures.

Cheers,
Rudi

> --
>
> ---
> You received this message because you are subscribed to the Google Groups "sage-matroid" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sage-matroid...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.

[Dima Pasechnik]

On 2013-08-27, Rudi Pendavingh <rudi.pe...@gmail.com> wrote:
> Hi Dima,
>
> I don't think we discuss the implementation in much detail in the documentation strings.
>
> In case you wonder, the algorithm we implemented is *not* the poly-time algorithm that follows from regular matroid decomposition. It is a pretty generic algorithm that also serves to test if a matroid is representable over other partial fields.
>
> To test if the matrix that represents a matroid M is TU, we consider each flat F of rank r(M)-2. Contracting such a flat in the represented matroid M, we get a represented minor N=M/F of rank 2. Say this N is represented by the matrix
> [ 1 0 1 .. 1 ..]
> [ 0 1 1 .. a_i..]
> If one of the a_i is not in {0,1}, then the representation of M is not TU. Conversely, if M passes this test for each F, then its representation is TU.

Looks puzzling for someone who never dealt with matroid minors before :)
Surely there must be some obvious linear algebra behind this, no?

>
> So it is an inherently exponential algorithm, but fairly optimized in its use of data structures.

How big is the rank it can realistically handle?
Would, say, 200 be out of question?

It might be I'll need something faster than this for my purposes. A
quick search brings up an implementation, using C++ and boost:
http://www.utdallas.edu/~klaus/TUtest/
of a simplified (and slower, at least in theory)
version of Truemper's algorithm.
Would it be worthwhile to create a Cython interface for it?
Or are there better codes to use?

Thanks,
Dima

[Rudi Pendavingh]

>> To test if the matrix that represents a matroid M is TU, we consider each flat F of rank r(M)-2. Contracting such a flat in the represented matroid M, we get a represented minor N=M/F of rank 2. Say this N is represented by the matrix
>> [ 1 0 1 .. 1 ..]
>> [ 0 1 1 .. a_i..]
>> If one of the a_i is not in {0,1}, then the representation of M is not TU. Conversely, if M passes this test for each F, then its representation is TU.
>
> Looks puzzling for someone who never dealt with matroid minors before :)
> Surely there must be some obvious linear algebra behind this, no?
>

There is. To naively test if A is totally unimodular, we need to calculate the determinants of all square submatrices A[X,Y]. If A' arises from A by pivoting in the xy-th entry, and x in X, y in Y, then A[X,Y]=A[X-x, Y-y]. Also, det A[X,Y] takes constant time to calculate if X and Y are of bounded size. So you can imagine that by cleverly pivoting around in A you will eventually have each subdeterminant 'close by' as the determinant of a small submatrix. The flats of ranks 2 in the above story is how we organize the pivots.

>>
>> So it is an inherently exponential algorithm, but fairly optimized in its use of data structures.
> How big is the rank it can realistically handle?
> Would, say, 200 be out of question?

Absolutely impossible. 20 already sounds hard.

>
> It might be I'll need something faster than this for my purposes. A
> quick search brings up an implementation, using C++ and boost:
> http://www.utdallas.edu/~klaus/TUtest/
> of a simplified (and slower, at least in theory)
> version of Truemper's algorithm.
> Would it be worthwhile to create a Cython interface for it?

A cython interface for any such code would be fantastic, Dima.

> Or are there better codes to use?

I do not know any other implementation than Truempers.
>

Cheers,
Rudi

[Dima Pasechnik]

On 2013-08-27, Rudi Pendavingh <rudi.pe...@gmail.com> wrote:
>
>

A mild problem is that it needs an older version of Boost C++ library
than the one in Sage. From Boost versions 1.51 and later
(more precisely, this commit into Boost
http://lists.boost.org/boost-commit/2012/03/40240.php)
broke the code, Boost graph library is rather esoteric
and I'd hate to fix it myself.
(the implementation uses very fancy C++ features I didn't even see
before, namely
http://www.boost.org/doc/libs/1_54_0/doc/html/variant.html)

Dima

[Volker Braun]

It just requires some trivial fixes, boost changed to long instead of int (really the TUtest library should use typedefs of boost instead of hardcoding int or long). Anyways:

--- src/graphicness.hpp~ 2012-10-12 13:20:30.000000000 +0100

+++ src/graphicness.hpp 2013-08-28 19:23:08.999432258 +0100

@@ -345,7 +345,7 @@

       typedef std::vector<vertex_t> vertex_vector_t;

       typedef std::vector<edge_t> edge_vector_t;

 

-      typedef boost::vec_adj_list_vertex_id_map<boost::no_property, unsigned int> IndexMap;

+      typedef boost::vec_adj_list_vertex_id_map<boost::no_property, unsigned long> IndexMap;

       IndexMap index_map = boost::get(boost::vertex_index, graph);

 

       /// Collect all 1-edges

@@ -474,7 +474,7 @@

         if (abort)

           continue;

 

-        boost::one_bit_color_map<boost::vec_adj_list_vertex_id_map<boost::no_property, unsigned int> > bipartition(

+        boost::one_bit_color_map<boost::vec_adj_list_vertex_id_map<boost::no_property, unsigned long> > bipartition(

             num_components, boost::get(boost::vertex_index, component_graph));

 

         if (!boost::is_bipartite(component_graph, boost::get(boost::vertex_index, component_graph), bipartition))

[Dima Pasechnik]

Hi Volker,

On 2013-08-28, Volker Braun <vbrau...@gmail.com> wrote:
> It just requires some trivial fixes, boost changed to long instead of int
> (really the TUtest library should use typedefs of boost instead of
> hardcoding int or long). Anyways:

it's surely a step in the right direction, but 'make check' fails with
similarly obscure error messages. That is, with gcc 4.7.
With gcc 4.2 things compile.

The culprit seems to be the following fragment in
src/matrix_transposed.hpp:

template <typename MatrixType>
inline void matrix_permute1(matrix_transposed <MatrixType>& matrix,
size_t index1, size_t index2)
{
matrix_permute2(matrix.data(), index1, index2);
}
[...]
template <typename MatrixType>
inline void matrix_permute2(matrix_transposed <MatrixType>& matrix,
size_t index1, size_t index2)
{
matrix_permute1(matrix.data(), index1, index2);
}


Any idea? I tried following an advice on porting to gcc 4.7,
forward-declaring matrix_permute1 and matrix_permute2, but got
nowhere.

Dima

[Volker Braun]

On Thursday, August 29, 2013 6:16:31 AM UTC+1, Dima Pasechnik wrote:

Any idea? I tried following an advice on porting to gcc 4.7, 

forward-declaring matrix_permute1 and matrix_permute2, but got
nowhere.

No, obviously this can't be solved by defining matrix_permut1 in terms of matrix_permute2 and defining matrix_permute2 in terms of matrix_permute1. There are functions with the same name that operate on matrix proxies and these are the ones that the compilation unit fails to see. Its another easy fix to declare them in the right order.

A fixed version is here:  https://github.com/vbraun/unimodularity-library

I do get segfaults if I compile it with -O2 and higher, but

./configure CXXFLAGS='-O1' && make && make check

works.

[Dima Pasechnik]

On 2013-08-29, Volker Braun <vbrau...@gmail.com> wrote:
> On Thursday, August 29, 2013 6:16:31 AM UTC+1, Dima Pasechnik wrote:
>
>> Any idea? I tried following an advice on porting to gcc 4.7,
>>
> forward-declaring matrix_permute1 and matrix_permute2, but got
>> nowhere.
>
>
> No, obviously this can't be solved by defining matrix_permut1 in terms of
> matrix_permute2 and defining matrix_permute2 in terms of matrix_permute1.
> There are functions with the same name that operate on matrix proxies and
> these are the ones that the compilation unit fails to see. Its another easy
> fix to declare them in the right order.

Same name?! If so, OMFG... It beats me.
It's been over 10 years since I touched C++, and
I better stay away from it as much as I can :-)

>
> A fixed version is here: https://github.com/vbraun/unimodularity-library

thanks!

[Volker Braun]

On Thursday, August 29, 2013 1:35:04 PM UTC+1, Dima Pasechnik wrote:

Same name?! If so, OMFG... It beats me.

Polymorphism is a feature of C++. You want to swap rows/columns of actual matrices (by copying the values around) and of matrix proxies (where you just change the map of proxied index -> actual index). So you use the same function name, and the implementation is chosen by the compiler depending on the type. Basically, its the C++ version of 

if isinstance(arg, TypeA):

    ...

elif isinstance(arg, TypeB):

   ...

[Dima Pasechnik]

On 2013-08-29, Volker Braun <vbrau...@gmail.com> wrote:

I'm perfectly aware of template functionality in C++, derived classes,
etc, but you wrote that matrix_permute1 and matrix_permute2
are the same function!

> > There are functions with the same name

Is this true? Boost/variant got be so scared that
I'd not be suprised if it was true...

Dima

[Dima Pasechnik]

I miscounted the number of different implementations of
matrix_permute1 and matrix_permute2 in the source. :-)
Sorry for noise...
>
> Dima
>
>