Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Linear implementation for Fulkerson-Chen-Anstee #2537

Merged
merged 17 commits into from Apr 1, 2024

Conversation

Tagl
Copy link
Contributor

@Tagl Tagl commented Mar 18, 2024

Closes #1452

The current implementation of the Fulkerson-Chen-Anstee simple directed graphicality check is quadratic.
This is my own approach to the linear time version of Fulkerson-Chen-Anstee.
I (with @szhorvat) plan to do a write-up for this approach.

After implementing the method I read section 5.1 in:
H. Kim, C. I. Del Genio, K. E. Bassler, and Z. Toroczkai, New J. Phys. 14, 023012 (2012). https://dx.doi.org/10.1088/1367-2630/14/2/023012
The methods follow the same general approach.
The key difference is they use recurrence relations to prove that the method works.
I leave here an explanation of how my method works and as in the other methods, the first step is sorting.

The theorem (screenshotted from Wikipedia):

image

The left hand side is easy to compute in linear time for each value of k.

On the right hand side we have two sums of outdegrees.
For the first sum we separate outdegrees based on whether they are $< k-1$ or $>= k-1$.
Similarly, for the second sum we separate outdegrees based on whether they are $< k$ or $>= k$.
In both cases we have a value which we can think of as the splitting point of a sum.
We only need to track the sum of outdegrees below the splitting point, because the value of k is incrementing from $1$ to $n$, and therefore they will not individually be impacted by future changes.

The outdegrees that are equal to or above the splitting point must be tracked, and we do so with an ordered multiset.
The ordered set for the first sum is initially empty and the one for the second sum initially contains all of the outdegrees.
As $k$ increases, either add/remove the outdegree at index k-1 to/from the sum or ordered set, depending on which side of the splitting point it lies.
We may also need to remove the smallest outdegrees from the ordered sets and move them over to the sums of outdegrees below the splitting point.
The contribution of the values in the ordered multiset to the right hand side is the value of $k-1$ or $k$, depending on which sum we are considering, multiplied by the size of the ordered multiset.

Using a balanced binary search tree for the ordered set means adding and removing from the set takes $\mathcal{O}(\log n)$ time.
Then the time complexity will match that of the sort operation.

We can push this a little further.
Since the values in the set are bounded by n, we can switch the quicksort out for a bucket sort which is $\mathcal{O}(n)$.
Additionally, it allows us to switch out a balanced binary search tree for an array of integers, where the index in to the array represents the degree and the value at that index represents how many times that degree appears in the ordered multiset.
All operations within the main loop are $\mathcal{O}(1)$ and therefore the time complexity of the algorithm is $\mathcal{O}(n)$ in total.
I have not done any benchmarking yet, but I have tested an equivalent C++ implementation on sequences of size $1$ to $2,000,000$ and there is a noticeable difference in running time between the linear and the linearithmic approaches.

@Tagl Tagl force-pushed the linear_fulkerson_chen_anstee branch from d4fcc86 to 7b46439 Compare March 18, 2024 16:16
Copy link

codecov bot commented Mar 18, 2024

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 84.38%. Comparing base (808c083) to head (bb55b15).
Report is 60 commits behind head on master.

❗ Current head bb55b15 differs from pull request most recent head 05637a2. Consider uploading reports for the commit 05637a2 to get more accurate results

Additional details and impacted files

Impacted file tree graph

@@           Coverage Diff           @@
##           master    #2537   +/-   ##
=======================================
  Coverage   84.37%   84.38%           
=======================================
  Files         376      376           
  Lines       61633    61691   +58     
  Branches    12060    12074   +14     
=======================================
+ Hits        52004    52059   +55     
- Misses       9629     9632    +3     
Files Coverage Δ
src/misc/graphicality.c 98.00% <100.00%> (+0.26%) ⬆️

Continue to review full report in Codecov by Sentry.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update af2d95e...05637a2. Read the comment docs.

Copy link
Member

@ntamas ntamas left a comment

Choose a reason for hiding this comment

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

A few comments, nothing conclusive yet. I've talked to @szhorvat and he will also want to chime in and review later.

src/misc/graphicality.c Outdated Show resolved Hide resolved
src/misc/graphicality.c Show resolved Hide resolved
@szhorvat
Copy link
Member

I edited the theorem into the top post so it's easier to follow your description, and what the LHS and RHS refer to.

@szhorvat
Copy link
Member

szhorvat commented Mar 31, 2024

I added some rudimentary benchmarks.

EDIT: First version of the benchmark was buggy, now fixed the results below.


Before this PR:

|> Benchmark file: /Users/szhorvat/Repos/igraph-main/igraph/tests/benchmarks/graphicality.c
| Undirected simple,       PA,     n = 50, 2000000x                                 0.22s   0.22s      0s
| Undirected loops,        PA,     n = 50, 2000000x                                0.596s  0.594s      0s
| Undirected multi,        PA,     n = 50, 2000000x                                 0.05s  0.049s      0s
| Undirected multi, loops, PA,     n = 50, 2000000x                                    0s      0s      0s
| Undirected simple,       G(n,p), n = 50, 2000000x                                0.182s   0.18s      0s
| Undirected loops,        G(n,p), n = 50, 2000000x                                0.775s  0.772s      0s
| Undirected multi,        G(n,p), n = 50, 2000000x                                 0.05s  0.049s      0s
| Undirected multi, loops, G(n,p), n = 50, 2000000x                                    0s      0s      0s
| Directed simple,       PA,     n = 50, 2000000x                                   1.26s   1.26s      0s
| Directed loops,        PA,     n = 50, 2000000x                                   1.15s   1.15s      0s
| Directed multi,        PA,     n = 50, 2000000x                                  0.062s  0.061s      0s
| Directed multi, loops, PA,     n = 50, 2000000x                                      0s      0s      0s
| Directed simple,       G(n,p), n = 50, 2000000x                                   1.71s   1.71s      0s
| Directed loops,        G(n,p), n = 50, 2000000x                                   1.51s   1.51s      0s
| Directed multi,        G(n,p), n = 50, 2000000x                                  0.062s  0.061s      0s
| Directed multi, loops, G(n,p), n = 50, 2000000x                                      0s      0s      0s

| Undirected simple,       PA,     n = 1000, 100000x                               0.128s  0.126s      0s
| Undirected loops,        PA,     n = 1000, 100000x                               0.539s  0.536s      0s
| Undirected multi,        PA,     n = 1000, 100000x                               0.058s  0.057s      0s
| Undirected multi, loops, PA,     n = 1000, 100000x                                   0s      0s      0s
| Undirected simple,       G(n,p), n = 1000, 100000x                               0.101s    0.1s      0s
| Undirected loops,        G(n,p), n = 1000, 100000x                               0.673s  0.671s      0s
| Undirected multi,        G(n,p), n = 1000, 100000x                               0.058s  0.057s      0s
| Undirected multi, loops, G(n,p), n = 1000, 100000x                                   0s      0s      0s
| Directed simple,       PA,     n = 1000, 100000x                                  1.89s   1.89s      0s
| Directed loops,        PA,     n = 1000, 100000x                                 0.822s  0.818s      0s
| Directed multi,        PA,     n = 1000, 100000x                                  0.06s  0.059s      0s
| Directed multi, loops, PA,     n = 1000, 100000x                                     0s      0s      0s
| Directed simple,       G(n,p), n = 1000, 100000x                                  3.56s   3.54s  0.001s
| Directed loops,        G(n,p), n = 1000, 100000x                                   1.9s   1.89s      0s
| Directed multi,        G(n,p), n = 1000, 100000x                                  0.06s   0.06s      0s
| Directed multi, loops, G(n,p), n = 1000, 100000x                                     0s      0s      0s

| Undirected simple,       PA,     n = 1000000, 100x                               0.151s  0.133s  0.016s
| Undirected loops,        PA,     n = 1000000, 100x                                1.01s  0.989s   0.02s
| Undirected multi,        PA,     n = 1000000, 100x                               0.058s  0.058s      0s
| Undirected multi, loops, PA,     n = 1000000, 100x                               0.059s  0.057s      0s
| Undirected simple,       G(n,p), n = 1000000, 100x                                0.12s  0.102s  0.017s
| Undirected loops,        G(n,p), n = 1000000, 100x                                1.32s   1.29s  0.022s
| Undirected multi,        G(n,p), n = 1000000, 100x                               0.058s  0.058s      0s
| Undirected multi, loops, G(n,p), n = 1000000, 100x                               0.058s  0.057s      0s
| Directed simple,       PA,     n = 1000000, 100x                                  19.8s   19.7s  0.027s
| Directed loops,        PA,     n = 1000000, 100x                                   1.4s   1.35s  0.042s
| Directed multi,        PA,     n = 1000000, 100x                                 0.061s  0.061s      0s
| Directed multi, loops, PA,     n = 1000000, 100x                                 0.059s  0.058s      0s
| Directed simple,       G(n,p), n = 1000000, 100x                                  7.56s    7.5s  0.024s
| Directed loops,        G(n,p), n = 1000000, 100x                                  2.85s   2.79s  0.044s
| Directed multi,        G(n,p), n = 1000000, 100x                                 0.061s   0.06s      0s
| Directed multi, loops, G(n,p), n = 1000000, 100x                                 0.059s  0.058s      0s

After this PR:

|> Benchmark file: /Users/szhorvat/Repos/igraph-main/igraph/tests/benchmarks/graphicality.c
| Undirected simple,       PA,     n = 50, 2000000x                                0.232s  0.231s      0s
| Undirected loops,        PA,     n = 50, 2000000x                                0.669s  0.666s      0s
| Undirected multi,        PA,     n = 50, 2000000x                                0.055s  0.055s      0s
| Undirected multi, loops, PA,     n = 50, 2000000x                                    0s      0s      0s
| Undirected simple,       G(n,p), n = 50, 2000000x                                0.183s  0.182s      0s
| Undirected loops,        G(n,p), n = 50, 2000000x                                0.858s  0.856s      0s
| Undirected multi,        G(n,p), n = 50, 2000000x                                0.054s  0.053s      0s
| Undirected multi, loops, G(n,p), n = 50, 2000000x                                    0s      0s      0s
| Directed simple,       PA,     n = 50, 2000000x                                    3.1s   3.08s  0.001s
| Directed loops,        PA,     n = 50, 2000000x                                   1.24s   1.23s      0s
| Directed multi,        PA,     n = 50, 2000000x                                  0.063s  0.062s      0s
| Directed multi, loops, PA,     n = 50, 2000000x                                      0s      0s      0s
| Directed simple,       G(n,p), n = 50, 2000000x                                   3.29s   3.27s  0.001s
| Directed loops,        G(n,p), n = 50, 2000000x                                   1.61s   1.61s      0s
| Directed multi,        G(n,p), n = 50, 2000000x                                  0.063s  0.063s      0s
| Directed multi, loops, G(n,p), n = 50, 2000000x                                      0s      0s      0s

| Undirected simple,       PA,     n = 1000, 100000x                               0.129s  0.127s      0s
| Undirected loops,        PA,     n = 1000, 100000x                                0.57s  0.566s      0s
| Undirected multi,        PA,     n = 1000, 100000x                               0.061s  0.061s      0s
| Undirected multi, loops, PA,     n = 1000, 100000x                                   0s      0s      0s
| Undirected simple,       G(n,p), n = 1000, 100000x                               0.107s  0.107s      0s
| Undirected loops,        G(n,p), n = 1000, 100000x                               0.707s  0.703s      0s
| Undirected multi,        G(n,p), n = 1000, 100000x                               0.058s  0.057s      0s
| Undirected multi, loops, G(n,p), n = 1000, 100000x                                   0s      0s      0s
| Directed simple,       PA,     n = 1000, 100000x                                   2.4s   2.35s  0.001s
| Directed loops,        PA,     n = 1000, 100000x                                 0.966s  0.965s      0s
| Directed multi,        PA,     n = 1000, 100000x                                 0.063s  0.063s      0s
| Directed multi, loops, PA,     n = 1000, 100000x                                     0s      0s      0s
| Directed simple,       G(n,p), n = 1000, 100000x                                  2.19s   2.17s      0s
| Directed loops,        G(n,p), n = 1000, 100000x                                  1.58s   1.58s      0s
| Directed multi,        G(n,p), n = 1000, 100000x                                 0.064s  0.063s      0s
| Directed multi, loops, G(n,p), n = 1000, 100000x                                     0s      0s      0s

| Undirected simple,       PA,     n = 1000000, 100x                               0.157s  0.135s  0.022s
| Undirected loops,        PA,     n = 1000000, 100x                                1.08s   1.05s   0.02s
| Undirected multi,        PA,     n = 1000000, 100x                               0.059s  0.058s      0s
| Undirected multi, loops, PA,     n = 1000000, 100x                               0.059s  0.057s      0s
| Undirected simple,       G(n,p), n = 1000000, 100x                               0.124s  0.104s  0.019s
| Undirected loops,        G(n,p), n = 1000000, 100x                                1.39s   1.35s  0.024s
| Undirected multi,        G(n,p), n = 1000000, 100x                               0.061s  0.061s      0s
| Undirected multi, loops, G(n,p), n = 1000000, 100x                               0.059s  0.058s      0s
| Directed simple,       PA,     n = 1000000, 100x                                  2.64s   2.23s  0.385s
| Directed loops,        PA,     n = 1000000, 100x                                  1.53s   1.46s  0.059s
| Directed multi,        PA,     n = 1000000, 100x                                 0.062s  0.062s      0s
| Directed multi, loops, PA,     n = 1000000, 100x                                 0.061s   0.06s      0s
| Directed simple,       G(n,p), n = 1000000, 100x                                  2.54s   2.15s  0.299s
| Directed loops,        G(n,p), n = 1000000, 100x                                  2.95s   2.89s  0.034s
| Directed multi,        G(n,p), n = 1000000, 100x                                 0.062s  0.062s      0s
| Directed multi, loops, G(n,p), n = 1000000, 100x                                  0.06s   0.06s      0s

@szhorvat
Copy link
Member

szhorvat commented Mar 31, 2024

It seems that at about 1000 vertices the old implementation is still doing better in some cases (PA) despite the bad complexity. This might be because there's a lot of memory allocation when using a list of vectors, and memory allocation is expensive. Although the system time (3rd column) doesn't increase significantly. Let's discuss next week what might be going on.

@szhorvat
Copy link
Member

szhorvat commented Apr 1, 2024

@Tagl I refactored the counting sort to use a single array instead of an array of arrays. Can you please have a look to make sure I didn't break anything? Efficiency is greatly improved now:

|> Benchmark file: /Users/szhorvat/Repos/igraph-main/igraph/tests/benchmarks/graphicality.c
| Undirected simple,       PA,     n = 50, 2000000x                                 0.24s  0.239s      0s
| Undirected loops,        PA,     n = 50, 2000000x                                0.597s  0.594s      0s
| Undirected multi,        PA,     n = 50, 2000000x                                 0.05s  0.049s      0s
| Undirected multi, loops, PA,     n = 50, 2000000x                                    0s      0s      0s
| Undirected simple,       G(n,p), n = 50, 2000000x                                0.178s  0.175s      0s
| Undirected loops,        G(n,p), n = 50, 2000000x                                0.775s  0.771s      0s
| Undirected multi,        G(n,p), n = 50, 2000000x                                 0.05s  0.049s      0s
| Undirected multi, loops, G(n,p), n = 50, 2000000x                                    0s      0s      0s
| Directed simple,       PA,     n = 50, 2000000x                                   1.26s   1.26s      0s
| Directed loops,        PA,     n = 50, 2000000x                                   1.14s   1.14s      0s
| Directed multi,        PA,     n = 50, 2000000x                                  0.062s  0.061s      0s
| Directed multi, loops, PA,     n = 50, 2000000x                                      0s      0s      0s
| Directed simple,       G(n,p), n = 50, 2000000x                                   1.11s    1.1s      0s
| Directed loops,        G(n,p), n = 50, 2000000x                                    1.5s    1.5s      0s
| Directed multi,        G(n,p), n = 50, 2000000x                                  0.062s  0.061s      0s
| Directed multi, loops, G(n,p), n = 50, 2000000x                                      0s      0s      0s

| Undirected simple,       PA,     n = 1000, 100000x                               0.128s  0.126s      0s
| Undirected loops,        PA,     n = 1000, 100000x                               0.532s  0.529s      0s
| Undirected multi,        PA,     n = 1000, 100000x                               0.058s  0.057s      0s
| Undirected multi, loops, PA,     n = 1000, 100000x                                   0s      0s      0s
| Undirected simple,       G(n,p), n = 1000, 100000x                               0.101s    0.1s      0s
| Undirected loops,        G(n,p), n = 1000, 100000x                               0.647s  0.645s      0s
| Undirected multi,        G(n,p), n = 1000, 100000x                               0.058s  0.057s      0s
| Undirected multi, loops, G(n,p), n = 1000, 100000x                                   0s      0s      0s
| Directed simple,       PA,     n = 1000, 100000x                                 0.812s  0.809s      0s
| Directed loops,        PA,     n = 1000, 100000x                                 0.816s  0.814s      0s
| Directed multi,        PA,     n = 1000, 100000x                                  0.06s  0.059s      0s
| Directed multi, loops, PA,     n = 1000, 100000x                                     0s      0s      0s
| Directed simple,       G(n,p), n = 1000, 100000x                                 0.583s  0.581s      0s
| Directed loops,        G(n,p), n = 1000, 100000x                                  1.28s   1.28s      0s
| Directed multi,        G(n,p), n = 1000, 100000x                                  0.06s   0.06s      0s
| Directed multi, loops, G(n,p), n = 1000, 100000x                                     0s      0s      0s

| Undirected simple,       PA,     n = 1000000, 100x                               0.152s  0.133s  0.018s
| Undirected loops,        PA,     n = 1000000, 100x                                1.02s  0.992s  0.021s
| Undirected multi,        PA,     n = 1000000, 100x                               0.058s  0.058s      0s
| Undirected multi, loops, PA,     n = 1000000, 100x                               0.059s  0.058s      0s
| Undirected simple,       G(n,p), n = 1000000, 100x                               0.121s  0.102s  0.018s
| Undirected loops,        G(n,p), n = 1000000, 100x                                1.32s   1.29s  0.023s
| Undirected multi,        G(n,p), n = 1000000, 100x                               0.058s  0.058s      0s
| Undirected multi, loops, G(n,p), n = 1000000, 100x                               0.058s  0.057s      0s
| Directed simple,       PA,     n = 1000000, 100x                                     1s  0.852s  0.138s
| Directed loops,        PA,     n = 1000000, 100x                                  1.41s   1.36s  0.044s
| Directed multi,        PA,     n = 1000000, 100x                                  0.06s   0.06s      0s
| Directed multi, loops, PA,     n = 1000000, 100x                                 0.058s  0.058s      0s
| Directed simple,       G(n,p), n = 1000000, 100x                                 0.791s  0.627s  0.163s
| Directed loops,        G(n,p), n = 1000000, 100x                                  2.84s   2.79s  0.045s
| Directed multi,        G(n,p), n = 1000000, 100x                                  0.06s   0.06s      0s
| Directed multi, loops, G(n,p), n = 1000000, 100x                                 0.059s  0.058s      0s

@szhorvat
Copy link
Member

szhorvat commented Apr 1, 2024

@ntamas I think we're almost good to go with this. Any other comments? I ended up not using the single exit point you suggested because it makes it difficult to de-allocate memory early when we don't need it.

@szhorvat
Copy link
Member

szhorvat commented Apr 1, 2024

In principle a further optimization is possible since we already compute the cumulative in-degree sums, which are the LHS, during the counting sort, except that this code does it in the wrong order (from smallest to largest instead of largest to smallest, which the lhs is supposed to be).

This should not block merging though, as it's not even clear if it improves performance. For the final version and explanation of the algorithm we should sort it out though.

@Tagl
Copy link
Contributor Author

Tagl commented Apr 1, 2024

I refactored the counting sort to use a single array instead of an array of arrays. Can you please have a look to make sure I didn't break anything?

Looks good to me. The benchmark numbers are also much more in line with the measurements of my C++ implementation.

@szhorvat szhorvat merged commit 2d455d8 into igraph:master Apr 1, 2024
18 of 23 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.

is_graphical: linear-time test for simple directed case
3 participants