Skip to content

Replace implementation of maximum_weighted_matching()#400

Merged
jeremy-murphy merged 18 commits intoboostorg:developfrom
jorisvr:matching_joris
Feb 10, 2025
Merged

Replace implementation of maximum_weighted_matching()#400
jeremy-murphy merged 18 commits intoboostorg:developfrom
jorisvr:matching_joris

Conversation

@jorisvr
Copy link
Contributor

@jorisvr jorisvr commented Dec 5, 2024

This is a re-implementation of maximum_weighted_matching, based on a paper by Zvi Galil.
The new code runs in time O(V^3).
A new set of testcases are also added.

Resolves #199 #223 #399

The code has been tested extensively on large random graphs, using LEMON as a reference.

Faster algorithms are known for this problem. I initially planned to implement the O(VElog(V)) algorithm by Galil, Micali and Gabow. However it needs a mutable heap with features that are not readily available in the BGL, and it needs a special kind of mergeable priority queue. While possible, I feel that the amount of code would be disproportionate. So I decided to fall back to a simpler O(V^3) algorithm, essentially the same algorithm that inspired the previous implementation.

Feedback is very welcome. I will already mention a few points that may draw criticism:

  • I kept brute_force_maximum_weighted_matching() unchanged. This function is not very useful in my opinion, but it was part of the public API and there is no need to change it.
  • The documented API of maximum_weighted_matching() is backwards compatible with the previous code. But I removed the class weighted_augmenting_path_finder, which was essentially an internal detail although it lived in the global boost namespace.
  • I do runtime checks on some aspects of the input graph (vertex indices and range of edge weights). I don't see this much in the BGL and I guess it may be frowned upon. The thing is, the code will trigger undefined behavior if these preconditions are violated, and I feel like I can't let that happen.
  • Once a matching has been computed, a separate (much faster) algorithm can verify that the matching is optimal. If the primary algorithm is correct, this verification will never fail. I enabled the verification step by default, even though it is redundant and never fails. The primary algorithm is tricky. I feel that the certainty provided by the verification step is worth more than the clock cycles it costs.

Without this check, the test program declares all tests passed if it
fails to open the input file.
- Hand-picked graphs to explore basic functionality.
- Graphs that are known to trigger bugs in the old
  implementation of maximum_weighted_matching().
- Random small graphs.
The new code runs in O(V^3).
It also solves a number of known issues.
@jorisvr jorisvr changed the title Add tests for maximum_weighted_matching() Replace implementation of maximum_weighted_matching() Jan 13, 2025
@jeremy-murphy jeremy-murphy self-assigned this Jan 13, 2025
@jeremy-murphy
Copy link
Contributor

Are the heap data structures in Boost.Heap not sufficient? They are mergeable and mutable.
However, even if they are sufficient, I would be happy to simply get a correct implementation first that fixes the bugs and get a fast implementation later.
What do you think?

@jorisvr
Copy link
Contributor Author

jorisvr commented Jan 14, 2025

Are the heap data structures in Boost.Heap not sufficient?

The mutable heap in Boost.Heap would work as the "plain" type of heap in the matching algo. But it looks like BGL currently does not use Boost.Heap and I don't know how you feel about adding that dependency.

I also need a concatenable heap which does not currently exist in Boost. The merge feature of Boost.Heap is not sufficient. I need to merge heaps in O(log(n)) time with the option to unmerge them later in O(log(n)). The typical way to implement this is with a custom balanced binary tree. It's not rocket science but it adds another 800 lines or so.

LEMON and LEDA implement the O(V E log(V)) matching algorithm. It is much faster than O(V^3) on certain classes of sparse graphs. The speedup on random sparse graphs is fairly modest in my experience. And it can be slower on dense graphs.

The new code is already an order of magnitude faster than the previous version for graphs with V > 200. My feeling is that the faster algorithm adds a lot of code in exchange for little benefit. But I'm up for the challenge. If you want the best matching algorithm in BGL, I will be happy to work on it.

@jeremy-murphy
Copy link
Contributor

Thanks for the explanation. There's no problem with adding Boost.Heap as a dependency, as Boost.Graph already depends on many other parts of Boost. Sounds like the efficient algorithm would require adding a new data structure to Boost.Heap to start with, which shouldn't be too difficult, although I'm aware that the maintainer is not all that active any more. Given that the new implementation is much faster anyway, let's defer the efficient algorithm to later. Ultimately it would be nice to have a top-level algorithm that uses a heuristic to pick the fastest algorithm but users are still free to call specific algorithms. (Best of both worlds.)
Ok, now I have to find time to look over this code. Honestly, it could take a couple of weeks, so please be patient.
Thanks for the work!

@jorisvr
Copy link
Contributor Author

jorisvr commented Jan 15, 2025

Sounds like the efficient algorithm would require adding a new data structure to Boost.Heap to start with

I think the concatenable queue may be so special-purpose that it could just stay in BGL, but generalizing it is definitely also a valid option.

Given that the new implementation is much faster anyway, let's defer the efficient algorithm to later.

Agreed. It occurs to me that the O(V E log(V)) algorithm also needs an edge_index property, which breaks backward compatibility. There may be ways to deal with this, but it seems like it will be a more difficult road than the current PR.

Ok, now I have to find time to look over this code.

I understand. There is no hurry from my side. Thanks for supporting this effort.

@jeremy-murphy
Copy link
Contributor

What's the "nearest ancestor" problem referred to in the documentation for the fast Gabow algorithm? Is that LCA or something else?

@jorisvr
Copy link
Contributor Author

jorisvr commented Jan 22, 2025

What's the "nearest ancestor" problem referred to in the documentation for the fast Gabow algorithm?

To be honest, I don't know. I kept this comment from the documentation by Yi Ji as I saw no reason to remove it.

I know about the existence of that fast algorithm. I tried to read the paper by Gabow but I can not make heads or tails of it. Mehlhorn and Schaefer made an offhand remark that this algorithm may be unpractical (https://dl.acm.org/doi/10.1145/944618.944622 page 7) but that was a long time ago. I'm not aware of any public available implementation.

@jeremy-murphy
Copy link
Contributor

Doesn't surprise me too much. Bender et. al. made a similar remark about the theoretically optimal algorithm for LCA: it just ain't worth it.

@jeremy-murphy
Copy link
Contributor

Copy link
Contributor

@jeremy-murphy jeremy-murphy left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I haven't even got to the code proper yet but here are a few requests and questions to start with.

Copy link
Contributor

@jeremy-murphy jeremy-murphy left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, just a trivial request, but since the library is C++14 now can you please use using instead of typedef and remove spaces from in-between nested template <<....>> brackets? Thanks. I will enable clang-format one day...
PS. I mean, a template should be written as
template <typename Foo>, etc.

Copy link
Contributor

@jeremy-murphy jeremy-murphy left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not quite finished...

typedef typename property_traits< VertexIndexMap >::value_type index_t;
typedef typename std::make_unsigned<index_t>::type unsigned_index_t;
auto nv = num_vertices(g);
std::vector<bool> got_vertex(nv);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Make sure this is really what you want, as opposed to a bitset.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I believe std::bitset<N> requires its size to be fixed at compile time. But I need an array which is sized at run time to match the number of vertices of the graph.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, sorry, I meant Boost's dynamic_bitset.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah ok. vector<bool> and dynamic_bitset both provide the functionality I need. I don't see a specific reason to prefer one or the other.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I forgot to address this one, whoops.
vector<bool> is generally avoided because of its unusual performance trying to satisfy the standard container interface (providing a reference to each element via iteration) whilst only using one bit per value, which requires the use of a proxy class, etc.
A dynamic bitset is just more honest about what it is and does.
It's not a big deal, might not need to be changed until someone makes some other changes in there.

template < typename Func >
static void for_vertices_in_blossom(const blossom_t* blossom, Func func)
{
const nontrivial_blossom_t* ntb = blossom->nontrivial();
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just use auto?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't understand.
Do you mean auto func instead of template <typename Func> ? That was not allowed before C++20.

Or do you mean auto ntb instead of const nontrivial_blossom_t* ntb ?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The latter.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok. I changed it to auto for these and a few similar verbose declarations.

Is that what you had in mind, or do you want to push further towards the almost-alway-auto style?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I personally use and recommend AAA style, but what you've done here is fine. I just felt those long typenames with qualifiers were not helping the readability.

Also remove unnecessary < > spaces around template arguments.
Copy link
Contributor

@jeremy-murphy jeremy-murphy left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Still going,..

blossom_label_t label;

/** True if this is an instance of nontrivial_blossom. */
const bool is_nontrivial_blossom;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No const member variables either, same reason as references, but I'll review again on a proper screen to decide if it's worth changing.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's a good general principle, so yeah, please change it.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done.

typedef typename property_traits< VertexIndexMap >::value_type index_t;
typedef typename std::make_unsigned<index_t>::type unsigned_index_t;
auto nv = num_vertices(g);
std::vector<bool> got_vertex(nv);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, sorry, I meant Boost's dynamic_bitset.

template < typename Func >
static void for_vertices_in_blossom(const blossom_t* blossom, Func func)
{
const nontrivial_blossom_t* ntb = blossom->nontrivial();
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The latter.

= sub_blossom->vertices();
for (vertex_vec_iter_t v = sub_vertices.begin();
v != sub_vertices.end(); ++v)
if ((! edge.has_value()) || (s < slack))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm personally not in favour of a space between a unary operator such as ! and its operand, same as for other unary operators like dereference *, negation -, etc.
I don't think BGL has a history of it, but maybe I'm wrong?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh, I probably came up with this. I don't see it anywhere else in BGL.
I removed the spaces now.

I really like the space, but this is not the time to start a discussion about code formatting.

@jeremy-murphy
Copy link
Contributor

There are some functions that take a blossom* but require that it is not null, in which case they should take a reference instead.

@jeremy-murphy
Copy link
Contributor

That's all from me, once the last comments are resolved I'll merge it in.
Overall I love the code, so thank you!

@jorisvr
Copy link
Contributor Author

jorisvr commented Feb 6, 2025

There are some functions that take a blossom* but require that it is not null, in which case they should take a reference instead.

I now changed these into blossom&.
This involved adding a number of * and & operators in the code. Because blossom* is still the native type in several data structures, I end up going back-and-forth between pointer and reference. In my opinion, it makes the code less clear. But I can live with it.

@jeremy-murphy
Copy link
Contributor

jeremy-murphy commented Feb 10, 2025

There are some functions that take a blossom* but require that it is not null, in which case they should take a reference instead.

I now changed these into blossom&. This involved adding a number of * and & operators in the code. Because blossom* is still the native type in several data structures, I end up going back-and-forth between pointer and reference. In my opinion, it makes the code less clear. But I can live with it.

I know what you mean, sometimes you just have to pay a price for doing things correctly, but it also might mean that those other places using blossom* could do with some refactoring.

PS. And I realize that sounds hypocritical because I asked you to change the member variables from reference to pointer. Sometimes C++ is just ugly.

@jeremy-murphy jeremy-murphy merged commit 167ac18 into boostorg:develop Feb 10, 2025
22 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

maximum_weighted_matching() segfaults with constant weights

2 participants