diff --git a/.gitignore b/.gitignore index 0958ff7..e151a35 100644 --- a/.gitignore +++ b/.gitignore @@ -5,4 +5,6 @@ /docs/build/ # /docs/Manifest.toml /docs/src/index.md -.vscode/ \ No newline at end of file +.vscode/ +/benchmark/Manifest.toml +/benchmark/.CondaPkg/ diff --git a/benchmark/CondaPkg.toml b/benchmark/CondaPkg.toml new file mode 100644 index 0000000..5a67aef --- /dev/null +++ b/benchmark/CondaPkg.toml @@ -0,0 +1,7 @@ +[deps] +python = "3.12.*" + +[pip.deps] +graspologic = "==3.4.4" +networkx = ">=3.0" +scipy = ">=1.10" diff --git a/benchmark/Project.toml b/benchmark/Project.toml new file mode 100644 index 0000000..4976739 --- /dev/null +++ b/benchmark/Project.toml @@ -0,0 +1,24 @@ +[deps] +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +CondaPkg = "992eb4ea-22a4-4c89-a5bb-47a3300528ab" +Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" +GraphsOptim = "e79ef3ae-79c8-4b97-b165-63f338db30c2" +HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" + +[sources] +GraphsOptim = {path = ".."} + +[compat] +CairoMakie = "0.15" +CondaPkg = "0.2" +Graphs = "1.15" +HiGHS = "1" +JuMP = "1" +PythonCall = "0.9" +julia = "1.11" diff --git a/benchmark/README.md b/benchmark/README.md new file mode 100644 index 0000000..02d7d71 --- /dev/null +++ b/benchmark/README.md @@ -0,0 +1,25 @@ +# Benchmarks + +Benchmarks of GraphsOptim.jl against other Julia and Python libraries (Graphs.jl, graspologic, NetworkX, SciPy), in terms of running time and solution quality. +The results are shown on the [Benchmarks](https://juliagraphs.org/GraphsOptim.jl/dev/benchmarks/) page of the documentation, which also describes the methodology. + +| File | Content | +| :--- | :--- | +| `run.jl` | entry point, runs the suites and draws the figures | +| `suites/*.jl` | one suite per family of problems: instances, methods and metrics | +| `common.jl` | shared tools: timing, calls to Python, CSV files | +| `plots.jl` | figures, saved in `docs/src/assets/benchmarks/` | +| `results/*.csv` | raw results | + +To update the results and figures, run from the root of the repository: + +```bash +julia --project=benchmark -e 'using Pkg; Pkg.instantiate()' +julia --project=benchmark benchmark/run.jl # all suites +julia --project=benchmark benchmark/run.jl cliques # a single suite +julia --project=benchmark benchmark/run.jl --plots-only # only redraw the figures +``` + +The Python libraries are installed in `benchmark/.CondaPkg` by [CondaPkg.jl](https://github.com/JuliaPy/CondaPkg.jl), and called through [PythonCall.jl](https://github.com/JuliaPy/PythonCall.jl). + +To add a suite, create `suites/.jl` with a function `run_()` writing its CSV files, add a `plot_()` function to `plots.jl`, register both in `SUITES` in `run.jl`, and add a section to `docs/src/benchmarks.md`. diff --git a/benchmark/common.jl b/benchmark/common.jl new file mode 100644 index 0000000..943bc46 --- /dev/null +++ b/benchmark/common.jl @@ -0,0 +1,158 @@ +# Tools shared by all benchmark suites. + +# All libraries use a single BLAS thread for a fair comparison. +# This must be set before Python and NumPy are loaded. +ENV["OPENBLAS_NUM_THREADS"] = "1" + +using Graphs +using GraphsOptim +using HiGHS +using JuMP +using LinearAlgebra +using PythonCall +using Random +using SparseArrays +using Statistics + +BLAS.set_num_threads(1) + +const RESULTS_DIR = joinpath(@__DIR__, "results") +const FIGURES_DIR = joinpath(dirname(@__DIR__), "docs", "src", "assets", "benchmarks") + +"Time limit (in seconds) given to the optimizers and to the exponential-time algorithms." +const TIME_LIMIT = 60.0 + +"HiGHS with the benchmark time limit." +const OPTIMIZER = optimizer_with_attributes(HiGHS.Optimizer, "time_limit" => TIME_LIMIT) + +## Timing + +""" + measure(f; min_time=1.0, max_reps=7) + +Run `f()` once (after compilation by a previous call) and repeat it until `min_time` seconds or `max_reps` runs are reached. +Return the median running time and the result of the last run, or `(NaN, nothing)` if `f` fails, e.g. because a time limit was hit. +""" +function measure(f; min_time=1.0, max_reps=7) + times = Float64[] + local result + try + while length(times) < max_reps && (isempty(times) || sum(times) < min_time) + t = @elapsed result = f() + push!(times, t) + end + catch e + e isa InterruptException && rethrow() + @warn "Run failed" exception = (e, catch_backtrace()) + return NaN, nothing + end + return median(times), result +end + +pyexec( + """ + import time + + def measure(f, args, min_time=1.0, max_reps=7, **kwargs): + times = [] + while len(times) < max_reps and (not times or sum(times) < min_time): + t = time.perf_counter() + result = f(*args, **kwargs) + times.append(time.perf_counter() - t) + times.sort() + n = len(times) + median = times[n // 2] if n % 2 else (times[n // 2 - 1] + times[n // 2]) / 2 + return median, result + """, + Main, +) + +""" + pymeasure(f, args...; kwargs...) + +Same as [`measure`](@ref) for the Python function `f`, timed on the Python side so that the conversion of the arguments is not counted. +""" +function pymeasure(f, args...; kwargs...) + t, result = pyeval("measure", Main)(f, pytuple(args); kwargs...) + return pyconvert(Float64, t), result +end + +## Conversions + +"Convert an undirected or directed `Graphs.jl` graph into a NetworkX graph with vertices `0:n-1`." +function to_networkx(g::AbstractGraph; vertex_weights=nothing) + nx = pyimport("networkx") + h = is_directed(g) ? nx.DiGraph() : nx.Graph() + h.add_nodes_from(0:(nv(g) - 1)) + h.add_edges_from([(src(e) - 1, dst(e) - 1) for e in edges(g)]) + if !isnothing(vertex_weights) + for v in vertices(g) + h.nodes[v - 1]["weight"] = vertex_weights[v] + end + end + return h +end + +## Results + +function write_csv(file, rows) + mkpath(RESULTS_DIR) + open(joinpath(RESULTS_DIR, file), "w") do io + println(io, join(keys(first(rows)), ",")) + for row in rows + println(io, join(values(row), ",")) + end + end + return nothing +end + +function read_csv(file) + lines = readlines(joinpath(RESULTS_DIR, file)) + header = Tuple(Symbol.(split(lines[1], ","))) + parsevalue(x) = something(tryparse(Int, x), tryparse(Float64, x), String(x)) + return [ + NamedTuple{header}(Tuple(parsevalue.(split(line, ",")))) for line in lines[2:end] + ] +end + +function write_environment() + commit = try + readchomp(`git -C $(@__DIR__) describe --always --dirty`) + catch + "unknown" + end + version(mod) = pyconvert(String, pyimport(mod).__version__) + open(joinpath(RESULTS_DIR, "environment.md"), "w") do io + println(io, "- CPU: ", Sys.cpu_info()[1].model, " (single BLAS thread)") + println(io, "- Julia ", VERSION, ", GraphsOptim.jl at commit `", commit, "`") + println( + io, + "- Graphs.jl ", + pkgversion(Graphs), + ", JuMP.jl ", + pkgversion(JuMP), + ", HiGHS.jl ", + pkgversion(HiGHS), + ) + println( + io, + "- Python ", + pyconvert(String, pyimport("platform").python_version()), + ", graspologic ", + version("graspologic"), + ", NetworkX ", + version("networkx"), + ", SciPy ", + version("scipy"), + ) + end + return nothing +end + +""" + too_slow(time) + +Whether a method should be skipped on larger instances, because it failed or took more than a tenth of the time limit. +Exponential-time algorithms written in Python cannot be interrupted, so they are stopped early. +""" +too_slow(time) = isnan(time) || time > TIME_LIMIT / 10 diff --git a/benchmark/plots.jl b/benchmark/plots.jl new file mode 100644 index 0000000..b66b3c4 --- /dev/null +++ b/benchmark/plots.jl @@ -0,0 +1,284 @@ +# Figures of the benchmark results, saved as SVG files for the documentation. + +using CairoMakie + +const SURFACE = colorant"#fcfcfb" +const TEXT_PRIMARY = colorant"#0b0b0b" +const TEXT_SECONDARY = colorant"#52514e" +const GRID = colorant"#e4e3df" + +# Categorical palette, in fixed order. The color of a series follows its library. +const BLUE, ORANGE, AQUA, YELLOW = "#2a78d6", "#eb6834", "#1baf7a", "#eda100" +const MAGENTA, GREEN, VIOLET, RED = "#e87ba4", "#008300", "#4a3aa7", "#e34948" + +"Color and line style of each series: heuristics share the color of their library, dashed." +const SERIES_STYLE = Dict( + "GraphsOptim.jl" => (BLUE, :solid), + "GraphsOptim.jl (native)" => (BLUE, :solid), + "GraphsOptim.jl (JuMP)" => (VIOLET, :solid), + "graspologic" => (ORANGE, :solid), + "NetworkX" => (ORANGE, :solid), + "NetworkX (approximation)" => (ORANGE, :dash), + "NetworkX (2-approximation)" => (ORANGE, :dash), + "Graphs.jl" => (AQUA, :solid), + "Graphs.jl (all maximal cliques)" => (AQUA, :solid), + "Graphs.jl (greedy)" => (AQUA, :dash), + "Graphs.jl (2-approximation)" => (AQUA, :dot), + "SciPy" => (YELLOW, :solid), + # graph edit distance formulations + "F1" => (BLUE, :solid), + "F1prime" => (AQUA, :solid), + "F1plus" => (YELLOW, :solid), + "F2" => (MAGENTA, :solid), + "F2minus" => (GREEN, :solid), + "F2plus" => (VIOLET, :solid), + "FORI" => (RED, :solid), +) + +const THEME = Theme(; + backgroundcolor=SURFACE, + textcolor=TEXT_PRIMARY, + fontsize=13, + Axis=(; + width=280, + height=220, + backgroundcolor=SURFACE, + xgridcolor=GRID, + ygridcolor=GRID, + xgridwidth=1, + ygridwidth=1, + xminorgridvisible=false, + yminorgridvisible=false, + topspinevisible=false, + rightspinevisible=false, + leftspinecolor=TEXT_SECONDARY, + bottomspinecolor=TEXT_SECONDARY, + xtickcolor=TEXT_SECONDARY, + ytickcolor=TEXT_SECONDARY, + xticklabelcolor=TEXT_SECONDARY, + yticklabelcolor=TEXT_SECONDARY, + titlesize=14, + titlefont=:bold, + ), + Legend=(; framevisible=false, labelcolor=TEXT_PRIMARY, backgroundcolor=SURFACE), +) + +""" + series!(ax, rows; x, y, group, order, halo) + +Plot one line per series of `rows` (grouped by `group`), with `x` and `y` given as functions of a row. +The series listed in `halo` are drawn first as wide translucent bands, so that the other series stay visible when they coincide. +""" +function series!(ax, rows; x, y, group, order=unique(group.(rows)), halo=()) + for name in sort(order; by=name -> !(name in halo)) + rs = sort(filter(r -> group(r) == name, rows); by=x) + xs, ys = Float64.(x.(rs)), Float64.(y.(rs)) + keep = .!isnan.(ys) + any(keep) || continue + color, linestyle = SERIES_STYLE[name] + if name in halo + lines!(ax, xs[keep], ys[keep]; label=name, color=(color, 0.35), linewidth=8) + else + scatterlines!( + ax, + xs[keep], + ys[keep]; + label=name, + color, + linestyle, + linewidth=2, + markersize=8, + strokecolor=SURFACE, + strokewidth=2, + ) + end + end + readable_ticks!(ax, rows; x, y) + return nothing +end + +function legend_below!(fig, ax; nbanks=1) + return Legend( + fig[end + 1, :], + ax; + orientation=:horizontal, + nbanks, + tellheight=true, + tellwidth=false, + ) +end + +function format_time(t) + t >= 1 && return "$(round(Int, t)) s" + t >= 1e-3 && return "$(round(Int, 1e3t)) ms" + return "$(round(Int, 1e6t)) µs" +end + +"Readable ticks on logarithmic axes: the actual sizes horizontally, decades of seconds vertically." +function readable_ticks!(ax, rows; x, y) + if ax.xscale[] === log10 + xs = sort(unique(x.(rows))) + ax.xticks = (Float64.(xs), string.(xs)) + end + if ax.yscale[] === log10 + ys = filter(!isnan, Float64.(y.(rows))) + decades = 10.0 .^ (floor(Int, log10(minimum(ys))):ceil(Int, log10(maximum(ys)))) + # with few decades, add intermediate ticks + ticks = if maximum(ys) / minimum(ys) < 300 + [d * m for d in decades for m in (1, 2, 5)] + else + decades + end + ax.yticks = (ticks, format_time.(ticks)) + end + return nothing +end + +seconds_axis(pos; kwargs...) = Axis(pos; yscale=log10, ylabel="time", kwargs...) + +function save_figure(name, fig) + resize_to_layout!(fig) + return save(joinpath(FIGURES_DIR, name), fig) +end + +function plot_graph_matching() + scaling = read_csv("graph_matching_scaling.csv") + quality = read_csv("graph_matching_quality.csv") + with_theme(THEME) do + fig = Figure() + for (k, algorithm) in enumerate(("FAQ", "GOAT")) + ax = seconds_axis( + fig[1, k]; title=algorithm, xscale=log10, xlabel="number of vertices n" + ) + series!( + ax, + filter(r -> r.algorithm == algorithm, scaling); + x=r -> r.n, + y=r -> r.time, + group=r -> r.library, + ) + k == 1 && legend_below!(fig, ax) + end + save_figure("graph_matching_time.svg", fig) + + fig = Figure() + for (k, algorithm) in enumerate(("FAQ", "GOAT")) + ax = Axis( + fig[1, k]; + title=algorithm, + xlabel="correlation ρ", + ylabel="accuracy", + limits=(nothing, (-0.05, 1.05)), + ) + rows = [ + (; + ρ, + library, + accuracy=mean( + r.accuracy for r in quality if + r.ρ == ρ && r.library == library && r.algorithm == algorithm + ), + ) for ρ in unique(r.ρ for r in quality) for + library in ("GraphsOptim.jl", "graspologic") + ] + series!( + ax, + rows; + x=r -> r.ρ, + y=r -> r.accuracy, + group=r -> r.library, + halo=("graspologic",), + ) + k == 1 && legend_below!(fig, ax) + end + save_figure("graph_matching_accuracy.svg", fig) + end +end + +function plot_edit_distance() + rows = read_csv("edit_distance.csv") + with_theme(THEME) do + fig = Figure() + ax = seconds_axis(fig[1, 1]; xlabel="number of vertices n", width=600) + order = [string.(nameof.(GED_FORMULATIONS)); "NetworkX"] + series!(ax, rows; x=r -> r.n, y=r -> r.time, group=r -> r.method, order) + ax.xticks = sort(unique(r.n for r in rows)) + legend_below!(fig, ax; nbanks=2) + save_figure("edit_distance_time.svg", fig) + end +end + +function plot_paths_flows() + rows = read_csv("paths_flows.csv") + with_theme(THEME) do + fig = Figure() + for (k, problem) in enumerate(("shortest path", "min-cost flow", "assignment")) + ax = seconds_axis( + fig[1, k]; title=problem, xscale=log10, xlabel="number of vertices n" + ) + series!( + ax, + filter(r -> r.problem == problem, rows); + x=r -> r.n, + y=r -> r.time, + group=r -> r.library, + ) + Legend(fig[2, k], ax; framevisible=false, tellwidth=false, valign=:top) + end + save_figure("paths_flows_time.svg", fig) + end +end + +function plot_cliques() + rows = read_csv("cliques.csv") + problems = ("maximum weight clique", "maximum independent set", "minimum vertex cover") + with_theme(THEME) do + fig = Figure() + for (k, problem) in enumerate(problems) + ax = seconds_axis( + fig[1, k]; title=problem, xscale=log10, xlabel="number of vertices n" + ) + series!( + ax, + filter(r -> r.problem == problem, rows); + x=r -> r.n, + y=r -> r.time, + group=r -> r.method, + ) + Legend(fig[2, k], ax; framevisible=false, tellwidth=false, valign=:top) + end + save_figure("cliques_time.svg", fig) + + # Quality of the heuristics, relative to the optimum found by the exact methods + fig = Figure() + for (k, problem) in enumerate(problems[2:3]) + ax = Axis( + fig[1, k]; + title=problem, + xscale=log10, + xlabel="number of vertices n", + ylabel="size / optimal size", + ) + prows = filter(r -> r.problem == problem && r.valid == "true", rows) + relative = [] + for n in unique(r.n for r in prows) + exact = [r.value for r in prows if r.n == n && r.exact == "true"] + isempty(exact) && continue + opt = problem == "maximum independent set" ? maximum(exact) : minimum(exact) + for r in prows + r.n == n && push!(relative, (; n, r.method, ratio=r.value / opt)) + end + end + series!( + ax, + relative; + x=r -> r.n, + y=r -> r.ratio, + group=r -> r.method, + halo=("NetworkX",), + ) + Legend(fig[2, k], ax; framevisible=false, tellwidth=false, valign=:top) + end + save_figure("cliques_quality.svg", fig) + end +end diff --git a/benchmark/results/cliques.csv b/benchmark/results/cliques.csv new file mode 100644 index 0000000..46b913d --- /dev/null +++ b/benchmark/results/cliques.csv @@ -0,0 +1,65 @@ +problem,n,method,exact,time,value,valid +maximum weight clique,20,GraphsOptim.jl,true,0.005229667,30,true +maximum weight clique,20,Graphs.jl (all maximal cliques),true,2.3208e-5,30,true +maximum weight clique,20,NetworkX,true,9.7625e-5,30,true +maximum weight clique,40,GraphsOptim.jl,true,0.008800583,50,true +maximum weight clique,40,Graphs.jl (all maximal cliques),true,0.00045525,50,true +maximum weight clique,40,NetworkX,true,0.000469417,50,true +maximum weight clique,60,GraphsOptim.jl,true,0.12414775,58,true +maximum weight clique,60,Graphs.jl (all maximal cliques),true,0.002720209,58,true +maximum weight clique,60,NetworkX,true,0.001590375,58,true +maximum weight clique,80,GraphsOptim.jl,true,0.472089917,62,true +maximum weight clique,80,Graphs.jl (all maximal cliques),true,0.01433425,62,true +maximum weight clique,80,NetworkX,true,0.003569958,62,true +maximum weight clique,100,GraphsOptim.jl,true,0.673361375,69,true +maximum weight clique,100,Graphs.jl (all maximal cliques),true,0.036809,69,true +maximum weight clique,100,NetworkX,true,0.006860375,69,true +maximum weight clique,150,GraphsOptim.jl,true,9.650090583,72,true +maximum weight clique,150,Graphs.jl (all maximal cliques),true,0.28616475,72,true +maximum weight clique,150,NetworkX,true,0.041234125,72,true +maximum weight clique,200,Graphs.jl (all maximal cliques),true,1.738086875,77,true +maximum weight clique,200,NetworkX,true,0.128025791,77,true +maximum independent set,25,GraphsOptim.jl,true,0.006057416,10,true +maximum independent set,25,NetworkX,true,0.000286667,10,true +maximum independent set,25,Graphs.jl (greedy),false,1.834e-6,9,true +maximum independent set,25,NetworkX (approximation),false,0.001815875,9,true +minimum vertex cover,25,GraphsOptim.jl,true,0.013261458,15,true +minimum vertex cover,25,Graphs.jl (greedy),false,1.75e-6,15,true +minimum vertex cover,25,Graphs.jl (2-approximation),false,8.33e-7,20,true +minimum vertex cover,25,NetworkX (2-approximation),false,1.2042e-5,19,true +maximum independent set,50,GraphsOptim.jl,true,0.022556167,21,true +maximum independent set,50,NetworkX,true,0.003917417,21,true +maximum independent set,50,Graphs.jl (greedy),false,3.833e-6,20,true +maximum independent set,50,NetworkX (approximation),false,0.012148959,19,true +minimum vertex cover,50,GraphsOptim.jl,true,0.026696042,29,true +minimum vertex cover,50,Graphs.jl (greedy),false,3.625e-6,31,true +minimum vertex cover,50,Graphs.jl (2-approximation),false,1.209e-6,42,true +minimum vertex cover,50,NetworkX (2-approximation),false,2.1667e-5,43,true +maximum independent set,100,GraphsOptim.jl,true,0.085403709,43,true +maximum independent set,100,NetworkX,true,0.156261833,43,true +maximum independent set,100,Graphs.jl (greedy),false,8.042e-6,42,true +maximum independent set,100,NetworkX (approximation),false,0.079771084,36,true +minimum vertex cover,100,GraphsOptim.jl,true,0.095529375,57,true +minimum vertex cover,100,Graphs.jl (greedy),false,6.792e-6,60,true +minimum vertex cover,100,Graphs.jl (2-approximation),false,2.0e-6,84,true +minimum vertex cover,100,NetworkX (2-approximation),false,3.9917e-5,77,true +maximum independent set,200,GraphsOptim.jl,true,1.72203775,87,true +maximum independent set,200,NetworkX,true,51.697063167,87,true +maximum independent set,200,Graphs.jl (greedy),false,1.5792e-5,84,true +maximum independent set,200,NetworkX (approximation),false,0.6465547705000001,71,true +minimum vertex cover,200,GraphsOptim.jl,true,1.40136925,113,true +minimum vertex cover,200,Graphs.jl (greedy),false,1.4875e-5,122,true +minimum vertex cover,200,Graphs.jl (2-approximation),false,3.833e-6,172,true +minimum vertex cover,200,NetworkX (2-approximation),false,7.6208e-5,156,true +maximum independent set,400,GraphsOptim.jl,true,29.717599917,176,true +maximum independent set,400,Graphs.jl (greedy),false,3.7542e-5,173,true +maximum independent set,400,NetworkX (approximation),false,5.02385325,142,true +minimum vertex cover,400,GraphsOptim.jl,true,26.432987583,224,true +minimum vertex cover,400,Graphs.jl (greedy),false,3.2375e-5,239,true +minimum vertex cover,400,Graphs.jl (2-approximation),false,7.875e-6,334,true +minimum vertex cover,400,NetworkX (2-approximation),false,0.000166875,319,true +maximum independent set,800,Graphs.jl (greedy),false,8.5875e-5,350,true +maximum independent set,800,NetworkX (approximation),false,44.833363959,294,true +minimum vertex cover,800,Graphs.jl (greedy),false,6.6834e-5,467,true +minimum vertex cover,800,Graphs.jl (2-approximation),false,1.4459e-5,658,true +minimum vertex cover,800,NetworkX (2-approximation),false,0.000316708,623,true diff --git a/benchmark/results/edit_distance.csv b/benchmark/results/edit_distance.csv new file mode 100644 index 0000000..1da7509 --- /dev/null +++ b/benchmark/results/edit_distance.csv @@ -0,0 +1,49 @@ +n,edges_G,edges_H,method,time,value +4,4,6,F1,0.00649825,2.0 +4,4,6,F1prime,0.005606125,2.0 +4,4,6,F1plus,0.006054917,2.0 +4,4,6,F2,0.004719584,2.0 +4,4,6,F2minus,0.040121958,2.0 +4,4,6,F2plus,0.005402125,2.0 +4,4,6,FORI,0.006566625,2.0 +4,4,6,NetworkX,0.000688791973516345,2.0 +6,3,3,F1,0.007026667,2.0 +6,3,3,F1prime,0.005403417,2.0 +6,3,3,F1plus,0.005600708,2.0 +6,3,3,F2,0.004718417,2.0 +6,3,3,F2minus,0.0144185,2.0 +6,3,3,F2plus,0.005456833,2.0 +6,3,3,FORI,0.006039417,2.0 +6,3,3,NetworkX,0.010165582993067801,2.0 +8,9,9,F1,0.317395416,2.0 +8,9,9,F1prime,0.067758583,2.000000000000002 +8,9,9,F1plus,0.013552916,1.9999999999999742 +8,9,9,F2,0.147045,1.9999999999982674 +8,9,9,F2minus,1.177653541,1.9999999999999805 +8,9,9,F2plus,0.036437959,2.000000000000016 +8,9,9,FORI,0.024888917,2.0000000000000053 +8,9,9,NetworkX,0.7807432090630755,2.0 +10,13,14,F1,2.255427833,2.9999999999998868 +10,13,14,F1prime,0.631977042,2.9999999999999583 +10,13,14,F1plus,0.036517667,3.000000000000247 +10,13,14,F2,0.341728458,3.0 +10,13,14,F2minus,12.634142,2.9999999999912204 +10,13,14,F2plus,0.037046459,2.9999999999999956 +10,13,14,FORI,0.0609545,3.0 +10,13,14,NetworkX,0.8282274169614539,3.0 +12,30,30,F1,NaN,NaN +12,30,30,F1prime,13.784473417,3.9999999999999973 +12,30,30,F1plus,0.213055416,3.9999999999991642 +12,30,30,F2,42.825320625,3.9999999999999636 +12,30,30,F2plus,0.460361541,4.000000000000417 +12,30,30,FORI,0.566252,4.0 +12,30,30,NetworkX,NaN,NaN +14,37,39,F1plus,0.485347,3.9999999999996083 +14,37,39,F2plus,1.306716458,3.999999999999851 +14,37,39,FORI,2.813300958,4.000000000003662 +16,47,46,F1plus,2.231851,4.999999999997288 +16,47,46,F2plus,4.56608075,5.000000000733629 +16,47,46,FORI,6.062009625,5.000000000000365 +18,57,59,F1plus,2.886359833,5.999999999999426 +18,57,59,F2plus,9.984754333,6.000000000000292 +20,73,73,F1plus,25.799572208,5.999999999994151 diff --git a/benchmark/results/environment.md b/benchmark/results/environment.md new file mode 100644 index 0000000..0347d34 --- /dev/null +++ b/benchmark/results/environment.md @@ -0,0 +1,4 @@ +- CPU: Apple M4 (single BLAS thread) +- Julia 1.11.7, GraphsOptim.jl at commit `7f42780-dirty` +- Graphs.jl 1.15.0, JuMP.jl 1.32.0, HiGHS.jl 1.26.0 +- Python 3.12.14, graspologic 3.4.4, NetworkX 3.7, SciPy 1.17.1 diff --git a/benchmark/results/graph_matching_quality.csv b/benchmark/results/graph_matching_quality.csv new file mode 100644 index 0000000..e825ec2 --- /dev/null +++ b/benchmark/results/graph_matching_quality.csv @@ -0,0 +1,241 @@ +ρ,seed,algorithm,library,time,accuracy,disagreements,disagreements_truth +0.5,1,FAQ,GraphsOptim.jl,0.023600458,0.006666666666666667,1336,1018 +0.5,1,FAQ,graspologic,0.034272082964889705,0.006666666666666667,1374,1018 +0.5,1,GOAT,GraphsOptim.jl,0.155964208,0.04666666666666667,1316,1018 +0.5,1,GOAT,graspologic,0.09579870896413922,0.04666666666666667,1324,1018 +0.5,2,FAQ,GraphsOptim.jl,0.019884334,0.0,1312,1010 +0.5,2,FAQ,graspologic,0.05440287501551211,0.013333333333333334,1306,1010 +0.5,2,GOAT,GraphsOptim.jl,0.131846,0.02,1312,1010 +0.5,2,GOAT,graspologic,0.12274516699835658,0.02,1316,1010 +0.5,3,FAQ,GraphsOptim.jl,0.018525333,0.0,1327,1013 +0.5,3,FAQ,graspologic,0.06961591693107039,0.006666666666666667,1319,1013 +0.5,3,GOAT,GraphsOptim.jl,0.167674875,0.013333333333333334,1289,1013 +0.5,3,GOAT,graspologic,0.10405933298170567,0.006666666666666667,1297,1013 +0.5,4,FAQ,GraphsOptim.jl,0.020940333,0.006666666666666667,1350,1034 +0.5,4,FAQ,graspologic,0.05503245803993195,0.006666666666666667,1364,1034 +0.5,4,GOAT,GraphsOptim.jl,0.115027916,0.006666666666666667,1334,1034 +0.5,4,GOAT,graspologic,0.1059612079989165,0.006666666666666667,1310,1034 +0.5,5,FAQ,GraphsOptim.jl,0.018287,0.013333333333333334,1259,965 +0.5,5,FAQ,graspologic,0.04908754210919142,0.006666666666666667,1265,965 +0.5,5,GOAT,GraphsOptim.jl,0.092924959,0.006666666666666667,1277,965 +0.5,5,GOAT,graspologic,0.10785570798907429,0.006666666666666667,1273,965 +0.5,6,FAQ,GraphsOptim.jl,0.015873625,0.0,1334,972 +0.5,6,FAQ,graspologic,0.0535116249229759,0.0,1334,972 +0.5,6,GOAT,GraphsOptim.jl,0.143383,0.006666666666666667,1330,972 +0.5,6,GOAT,graspologic,0.10572762507945299,0.013333333333333334,1318,972 +0.5,7,FAQ,GraphsOptim.jl,0.018472,0.0,1341,999 +0.5,7,FAQ,graspologic,0.052528750034980476,0.013333333333333334,1357,999 +0.5,7,GOAT,GraphsOptim.jl,0.145492709,0.0,1309,999 +0.5,7,GOAT,graspologic,0.10518287506420165,0.0,1319,999 +0.5,8,FAQ,GraphsOptim.jl,0.018636709,0.006666666666666667,1316,968 +0.5,8,FAQ,graspologic,0.0469584169331938,0.006666666666666667,1342,968 +0.5,8,GOAT,GraphsOptim.jl,0.148809834,0.013333333333333334,1308,968 +0.5,8,GOAT,graspologic,0.11038858303800225,0.013333333333333334,1306,968 +0.5,9,FAQ,GraphsOptim.jl,0.021637834,0.006666666666666667,1401,1061 +0.5,9,FAQ,graspologic,0.04853449994698167,0.006666666666666667,1413,1061 +0.5,9,GOAT,GraphsOptim.jl,0.14964475,0.013333333333333334,1387,1061 +0.5,9,GOAT,graspologic,0.10291779099497944,0.013333333333333334,1391,1061 +0.5,10,FAQ,GraphsOptim.jl,0.0197745,0.006666666666666667,1335,973 +0.5,10,FAQ,graspologic,0.04421933297999203,0.013333333333333334,1327,973 +0.5,10,GOAT,GraphsOptim.jl,0.147829208,0.013333333333333334,1315,973 +0.5,10,GOAT,graspologic,0.10804879199713469,0.006666666666666667,1317,973 +0.6,1,FAQ,GraphsOptim.jl,0.020412084,0.02,1347,783 +0.6,1,FAQ,graspologic,0.05402262497227639,0.006666666666666667,1349,783 +0.6,1,GOAT,GraphsOptim.jl,0.108847125,0.03333333333333333,1329,783 +0.6,1,GOAT,graspologic,0.0667668340029195,0.04,1329,783 +0.6,2,FAQ,GraphsOptim.jl,0.019326458,0.006666666666666667,1291,815 +0.6,2,FAQ,graspologic,0.05330941604916006,0.0,1291,815 +0.6,2,GOAT,GraphsOptim.jl,0.148613583,0.013333333333333334,1291,815 +0.6,2,GOAT,graspologic,0.0984279999975115,0.006666666666666667,1293,815 +0.6,3,FAQ,GraphsOptim.jl,0.020487958,0.0,1302,814 +0.6,3,FAQ,graspologic,0.050627042073756456,0.0,1288,814 +0.6,3,GOAT,GraphsOptim.jl,0.12249525,0.02666666666666667,1306,814 +0.6,3,GOAT,graspologic,0.08716883300803602,0.02,1310,814 +0.6,4,FAQ,GraphsOptim.jl,0.013775583,0.0,1370,818 +0.6,4,FAQ,graspologic,0.053754292079247534,0.0,1362,818 +0.6,4,GOAT,GraphsOptim.jl,0.139451083,0.013333333333333334,1316,818 +0.6,4,GOAT,graspologic,0.10492483293637633,0.013333333333333334,1316,818 +0.6,5,FAQ,GraphsOptim.jl,0.019257125,0.006666666666666667,1280,788 +0.6,5,FAQ,graspologic,0.05454729206394404,0.0,1312,788 +0.6,5,GOAT,GraphsOptim.jl,0.140568,0.013333333333333334,1276,788 +0.6,5,GOAT,graspologic,0.10011162492446601,0.013333333333333334,1274,788 +0.6,6,FAQ,GraphsOptim.jl,0.019996958,0.0,1331,789 +0.6,6,FAQ,graspologic,0.053916083998046815,0.02,1323,789 +0.6,6,GOAT,GraphsOptim.jl,0.165545916,0.02,1315,789 +0.6,6,GOAT,graspologic,0.11484583304263651,0.013333333333333334,1319,789 +0.6,7,FAQ,GraphsOptim.jl,0.019993833,0.0,1351,817 +0.6,7,FAQ,graspologic,0.050586082972586155,0.0,1325,817 +0.6,7,GOAT,GraphsOptim.jl,0.070942375,0.03333333333333333,1337,817 +0.6,7,GOAT,graspologic,0.10707383300177753,0.04,1337,817 +0.6,8,FAQ,GraphsOptim.jl,0.014196834,0.006666666666666667,1326,776 +0.6,8,FAQ,graspologic,0.037882125005126,0.013333333333333334,1332,776 +0.6,8,GOAT,GraphsOptim.jl,0.138755167,0.02,1294,776 +0.6,8,GOAT,graspologic,0.10512850002851337,0.02,1312,776 +0.6,9,FAQ,GraphsOptim.jl,0.020362125,0.02,1386,844 +0.6,9,FAQ,graspologic,0.05466379097197205,0.006666666666666667,1416,844 +0.6,9,GOAT,GraphsOptim.jl,0.149949458,0.013333333333333334,1394,844 +0.6,9,GOAT,graspologic,0.10825708403717726,0.02,1392,844 +0.6,10,FAQ,GraphsOptim.jl,0.018181334,0.006666666666666667,1354,768 +0.6,10,FAQ,graspologic,0.05406070896424353,0.013333333333333334,1324,768 +0.6,10,GOAT,GraphsOptim.jl,0.155060458,0.013333333333333334,1320,768 +0.6,10,GOAT,graspologic,0.06981133308727294,0.013333333333333334,1338,768 +0.7,1,FAQ,GraphsOptim.jl,0.019129083,0.013333333333333334,1355,597 +0.7,1,FAQ,graspologic,0.044289750047028065,0.006666666666666667,1363,597 +0.7,1,GOAT,GraphsOptim.jl,0.155278625,0.07333333333333333,1323,597 +0.7,1,GOAT,graspologic,0.08786391606554389,0.06666666666666667,1325,597 +0.7,2,FAQ,GraphsOptim.jl,0.019084792,0.03333333333333333,1280,584 +0.7,2,FAQ,graspologic,0.04509699996560812,0.02666666666666667,1318,584 +0.7,2,GOAT,GraphsOptim.jl,0.133717541,0.05333333333333334,1266,584 +0.7,2,GOAT,graspologic,0.0994867499684915,0.05333333333333334,1268,584 +0.7,3,FAQ,GraphsOptim.jl,0.015271833,0.02666666666666667,1304,612 +0.7,3,FAQ,graspologic,0.05245029192883521,0.013333333333333334,1298,612 +0.7,3,GOAT,GraphsOptim.jl,0.159401667,0.05333333333333334,1252,612 +0.7,3,GOAT,graspologic,0.08162037492729723,0.05333333333333334,1278,612 +0.7,4,FAQ,GraphsOptim.jl,0.0168855,0.02666666666666667,1365,637 +0.7,4,FAQ,graspologic,0.05350270797498524,0.02,1333,637 +0.7,4,GOAT,GraphsOptim.jl,0.114774542,0.03333333333333333,1347,637 +0.7,4,GOAT,graspologic,0.11332291702274233,0.03333333333333333,1345,637 +0.7,5,FAQ,GraphsOptim.jl,0.013903625,0.013333333333333334,1289,599 +0.7,5,FAQ,graspologic,0.05044474999886006,0.006666666666666667,1297,599 +0.7,5,GOAT,GraphsOptim.jl,0.141112625,0.0,1251,599 +0.7,5,GOAT,graspologic,0.10269887500908226,0.006666666666666667,1255,599 +0.7,6,FAQ,GraphsOptim.jl,0.020833875,0.04,1350,612 +0.7,6,FAQ,graspologic,0.05417491705156863,0.02,1344,612 +0.7,6,GOAT,GraphsOptim.jl,0.140240542,0.04666666666666667,1328,612 +0.7,6,GOAT,graspologic,0.10952845809515566,0.04666666666666667,1338,612 +0.7,7,FAQ,GraphsOptim.jl,0.013465125,0.013333333333333334,1335,601 +0.7,7,FAQ,graspologic,0.053821791894733906,0.03333333333333333,1325,601 +0.7,7,GOAT,GraphsOptim.jl,0.142482375,0.02666666666666667,1311,601 +0.7,7,GOAT,graspologic,0.11650583299342543,0.02,1307,601 +0.7,8,FAQ,GraphsOptim.jl,0.019887125,0.03333333333333333,1323,589 +0.7,8,FAQ,graspologic,0.05354908399749547,0.04,1303,589 +0.7,8,GOAT,GraphsOptim.jl,0.126728291,0.02666666666666667,1301,589 +0.7,8,GOAT,graspologic,0.11809800006449223,0.03333333333333333,1297,589 +0.7,9,FAQ,GraphsOptim.jl,0.020964708,0.006666666666666667,1398,614 +0.7,9,FAQ,graspologic,0.04656183300539851,0.0,1374,614 +0.7,9,GOAT,GraphsOptim.jl,0.103934666,0.06666666666666667,1372,614 +0.7,9,GOAT,graspologic,0.08860929193906486,0.06666666666666667,1372,614 +0.7,10,FAQ,GraphsOptim.jl,0.016800041,0.013333333333333334,1337,559 +0.7,10,FAQ,graspologic,0.040899959043599665,0.0,1333,559 +0.7,10,GOAT,GraphsOptim.jl,0.120704417,0.12,1313,559 +0.7,10,GOAT,graspologic,0.1251678749686107,0.1,1305,559 +0.8,1,FAQ,GraphsOptim.jl,0.0189135,0.02,1355,409 +0.8,1,FAQ,graspologic,0.064797209110111,0.06,1347,409 +0.8,1,GOAT,GraphsOptim.jl,0.01262575,1.0,409,409 +0.8,1,GOAT,graspologic,0.019436624948866665,1.0,409,409 +0.8,2,FAQ,GraphsOptim.jl,0.019604834,0.006666666666666667,1301,405 +0.8,2,FAQ,graspologic,0.0786475830245763,0.0,1293,405 +0.8,2,GOAT,GraphsOptim.jl,0.153747958,0.08666666666666667,1269,405 +0.8,2,GOAT,graspologic,0.11982379201799631,0.06666666666666667,1279,405 +0.8,3,FAQ,GraphsOptim.jl,0.02052425,0.013333333333333334,1304,402 +0.8,3,FAQ,graspologic,0.055299917003139853,0.0,1300,402 +0.8,3,GOAT,GraphsOptim.jl,0.012992458,1.0,402,402 +0.8,3,GOAT,graspologic,0.018190625007264316,1.0,402,402 +0.8,4,FAQ,GraphsOptim.jl,0.020789,0.02,1375,407 +0.8,4,FAQ,graspologic,0.04862158407922834,0.013333333333333334,1369,407 +0.8,4,GOAT,GraphsOptim.jl,0.013781458,1.0,407,407 +0.8,4,GOAT,graspologic,0.01807087496854365,1.0,407,407 +0.8,5,FAQ,GraphsOptim.jl,0.019194333,0.006666666666666667,1246,398 +0.8,5,FAQ,graspologic,0.05358845798764378,0.006666666666666667,1264,398 +0.8,5,GOAT,GraphsOptim.jl,0.033066625,1.0,398,398 +0.8,5,GOAT,graspologic,0.03255337500013411,1.0,398,398 +0.8,6,FAQ,GraphsOptim.jl,0.017287708,0.02,1343,419 +0.8,6,FAQ,graspologic,0.05537666694726795,0.006666666666666667,1357,419 +0.8,6,GOAT,GraphsOptim.jl,0.09639625,1.0,419,419 +0.8,6,GOAT,graspologic,0.06870416703168303,1.0,419,419 +0.8,7,FAQ,GraphsOptim.jl,0.01853775,0.006666666666666667,1308,398 +0.8,7,FAQ,graspologic,0.054290708038024604,0.013333333333333334,1320,398 +0.8,7,GOAT,GraphsOptim.jl,0.328810667,0.02,1312,398 +0.8,7,GOAT,graspologic,0.07690866698976606,0.02,1318,398 +0.8,8,FAQ,GraphsOptim.jl,0.01931975,0.006666666666666667,1327,377 +0.8,8,FAQ,graspologic,0.05557033303193748,0.006666666666666667,1329,377 +0.8,8,GOAT,GraphsOptim.jl,0.043273708,1.0,377,377 +0.8,8,GOAT,graspologic,0.043818249949254096,1.0,377,377 +0.8,9,FAQ,GraphsOptim.jl,0.020533458,0.013333333333333334,1394,424 +0.8,9,FAQ,graspologic,0.05285812495276332,0.03333333333333333,1416,424 +0.8,9,GOAT,GraphsOptim.jl,0.019097417,1.0,424,424 +0.8,9,GOAT,graspologic,0.021103874896652997,1.0,424,424 +0.8,10,FAQ,GraphsOptim.jl,0.019488917,0.006666666666666667,1346,376 +0.8,10,FAQ,graspologic,0.046237208996899426,0.02666666666666667,1324,376 +0.8,10,GOAT,GraphsOptim.jl,0.012546041,1.0,376,376 +0.8,10,GOAT,graspologic,0.01775174995418638,1.0,376,376 +0.9,1,FAQ,GraphsOptim.jl,0.016417792,0.013333333333333334,1361,207 +0.9,1,FAQ,graspologic,0.05386904196348041,0.04666666666666667,1375,207 +0.9,1,GOAT,GraphsOptim.jl,0.00852725,1.0,207,207 +0.9,1,GOAT,graspologic,0.011388708953745663,1.0,207,207 +0.9,2,FAQ,GraphsOptim.jl,0.014828209,0.006666666666666667,1302,210 +0.9,2,FAQ,graspologic,0.05359133402816951,0.04,1278,210 +0.9,2,GOAT,GraphsOptim.jl,0.009664416,1.0,210,210 +0.9,2,GOAT,graspologic,0.020490708062425256,1.0,210,210 +0.9,3,FAQ,GraphsOptim.jl,0.019854333,0.04,1263,207 +0.9,3,FAQ,graspologic,0.053879333077929914,0.013333333333333334,1299,207 +0.9,3,GOAT,GraphsOptim.jl,0.00580325,1.0,207,207 +0.9,3,GOAT,graspologic,0.015938582946546376,1.0,207,207 +0.9,4,FAQ,GraphsOptim.jl,0.021240041,0.013333333333333334,1345,211 +0.9,4,FAQ,graspologic,0.04407504200935364,0.05333333333333334,1371,211 +0.9,4,GOAT,GraphsOptim.jl,0.0050395,1.0,211,211 +0.9,4,GOAT,graspologic,0.011586250038817525,1.0,211,211 +0.9,5,FAQ,GraphsOptim.jl,0.0190465,0.02,1285,211 +0.9,5,FAQ,graspologic,0.03657662507612258,0.04,1279,211 +0.9,5,GOAT,GraphsOptim.jl,0.009017792,1.0,211,211 +0.9,5,GOAT,graspologic,0.01444699999410659,1.0,211,211 +0.9,6,FAQ,GraphsOptim.jl,0.015412,0.006666666666666667,1347,205 +0.9,6,FAQ,graspologic,0.045959874987602234,0.11333333333333333,1357,205 +0.9,6,GOAT,GraphsOptim.jl,0.004311125,1.0,205,205 +0.9,6,GOAT,graspologic,0.010526500060223043,1.0,205,205 +0.9,7,FAQ,GraphsOptim.jl,0.014786959,0.06666666666666667,1325,193 +0.9,7,FAQ,graspologic,0.05405679205432534,0.006666666666666667,1327,193 +0.9,7,GOAT,GraphsOptim.jl,0.009273583,1.0,193,193 +0.9,7,GOAT,graspologic,0.011516792001202703,1.0,193,193 +0.9,8,FAQ,GraphsOptim.jl,0.015503917,0.013333333333333334,1346,202 +0.9,8,FAQ,graspologic,0.04232966701965779,0.03333333333333333,1332,202 +0.9,8,GOAT,GraphsOptim.jl,0.009044459,1.0,202,202 +0.9,8,GOAT,graspologic,0.013600083999335766,1.0,202,202 +0.9,9,FAQ,GraphsOptim.jl,0.019668333,0.02,1403,215 +0.9,9,FAQ,graspologic,0.03754399996250868,0.02666666666666667,1393,215 +0.9,9,GOAT,GraphsOptim.jl,0.004002458,1.0,215,215 +0.9,9,GOAT,graspologic,0.011176749947480857,1.0,215,215 +0.9,10,FAQ,GraphsOptim.jl,0.01862075,0.02666666666666667,1324,204 +0.9,10,FAQ,graspologic,0.053759792004711926,0.04,1326,204 +0.9,10,GOAT,GraphsOptim.jl,0.007808541,1.0,204,204 +0.9,10,GOAT,graspologic,0.014430625014938414,1.0,204,204 +1.0,1,FAQ,GraphsOptim.jl,0.0066125,1.0,0,0 +1.0,1,FAQ,graspologic,0.014111250056885183,1.0,0,0 +1.0,1,GOAT,GraphsOptim.jl,0.005422709,1.0,0,0 +1.0,1,GOAT,graspologic,0.010709041962400079,1.0,0,0 +1.0,2,FAQ,GraphsOptim.jl,0.006954125,1.0,0,0 +1.0,2,FAQ,graspologic,0.019300333107821643,1.0,0,0 +1.0,2,GOAT,GraphsOptim.jl,0.007753458,1.0,0,0 +1.0,2,GOAT,graspologic,0.00889104197267443,1.0,0,0 +1.0,3,FAQ,GraphsOptim.jl,0.0065495,1.0,0,0 +1.0,3,FAQ,graspologic,0.01779554202221334,1.0,0,0 +1.0,3,GOAT,GraphsOptim.jl,0.006789916,1.0,0,0 +1.0,3,GOAT,graspologic,0.008640624932013452,1.0,0,0 +1.0,4,FAQ,GraphsOptim.jl,0.00764025,1.0,0,0 +1.0,4,FAQ,graspologic,0.014243457932025194,1.0,0,0 +1.0,4,GOAT,GraphsOptim.jl,0.0057815,1.0,0,0 +1.0,4,GOAT,graspologic,0.00863358296919614,1.0,0,0 +1.0,5,FAQ,GraphsOptim.jl,0.00491775,1.0,0,0 +1.0,5,FAQ,graspologic,0.026598916039802134,1.0,0,0 +1.0,5,GOAT,GraphsOptim.jl,0.007352959,1.0,0,0 +1.0,5,GOAT,graspologic,0.009037666954100132,1.0,0,0 +1.0,6,FAQ,GraphsOptim.jl,0.005155708,1.0,0,0 +1.0,6,FAQ,graspologic,0.014339791028760374,1.0,0,0 +1.0,6,GOAT,GraphsOptim.jl,0.005464708,1.0,0,0 +1.0,6,GOAT,graspologic,0.00872949999757111,1.0,0,0 +1.0,7,FAQ,GraphsOptim.jl,0.004913542,1.0,0,0 +1.0,7,FAQ,graspologic,0.013974000001326203,1.0,0,0 +1.0,7,GOAT,GraphsOptim.jl,0.007954125,1.0,0,0 +1.0,7,GOAT,graspologic,0.009695707936771214,1.0,0,0 +1.0,8,FAQ,GraphsOptim.jl,0.007371959,1.0,0,0 +1.0,8,FAQ,graspologic,0.012473082984797657,1.0,0,0 +1.0,8,GOAT,GraphsOptim.jl,0.007606709,1.0,0,0 +1.0,8,GOAT,graspologic,0.009102415991947055,1.0,0,0 +1.0,9,FAQ,GraphsOptim.jl,0.005056042,1.0,0,0 +1.0,9,FAQ,graspologic,0.018034124979749322,1.0,0,0 +1.0,9,GOAT,GraphsOptim.jl,0.006309208,1.0,0,0 +1.0,9,GOAT,graspologic,0.008603791007772088,1.0,0,0 +1.0,10,FAQ,GraphsOptim.jl,0.004493459,1.0,0,0 +1.0,10,FAQ,graspologic,0.01620624994393438,1.0,0,0 +1.0,10,GOAT,GraphsOptim.jl,0.00604625,1.0,0,0 +1.0,10,GOAT,graspologic,0.011293792049400508,1.0,0,0 diff --git a/benchmark/results/graph_matching_scaling.csv b/benchmark/results/graph_matching_scaling.csv new file mode 100644 index 0000000..5940f1d --- /dev/null +++ b/benchmark/results/graph_matching_scaling.csv @@ -0,0 +1,21 @@ +n,algorithm,library,time,accuracy,disagreements,disagreements_truth +50,FAQ,GraphsOptim.jl,0.001085834,0.44,75,21 +50,FAQ,graspologic,0.006419709068723023,1.0,21,21 +50,GOAT,GraphsOptim.jl,0.000756083,1.0,21,21 +50,GOAT,graspologic,0.004522875067777932,1.0,21,21 +100,FAQ,GraphsOptim.jl,0.006193167,0.05,541,105 +100,FAQ,graspologic,0.0268864999525249,0.04,547,105 +100,GOAT,GraphsOptim.jl,0.002770875,1.0,105,105 +100,GOAT,graspologic,0.008391584036871791,1.0,105,105 +200,FAQ,GraphsOptim.jl,0.040110541,0.01,2625,375 +200,FAQ,graspologic,0.1100711670005694,0.005,2653,375 +200,GOAT,GraphsOptim.jl,0.016171334,1.0,375,375 +200,GOAT,graspologic,0.019382541999220848,1.0,375,375 +500,FAQ,GraphsOptim.jl,0.355570125,0.004,17988,2282 +500,FAQ,graspologic,1.3242363330209628,0.01,18048,2282 +500,GOAT,GraphsOptim.jl,0.081232084,1.0,2282,2282 +500,GOAT,graspologic,0.24870108300819993,1.0,2282,2282 +1000,FAQ,GraphsOptim.jl,2.216433291,0.003,77328,9036 +1000,FAQ,graspologic,10.040262583992444,0.003,77474,9036 +1000,GOAT,GraphsOptim.jl,0.2503195835,0.055,88922,9036 +1000,GOAT,graspologic,0.6955229794839397,0.055,88922,9036 diff --git a/benchmark/results/paths_flows.csv b/benchmark/results/paths_flows.csv new file mode 100644 index 0000000..3c2ccef --- /dev/null +++ b/benchmark/results/paths_flows.csv @@ -0,0 +1,42 @@ +problem,n,library,time,value +shortest path,100,GraphsOptim.jl,0.027957459,94 +shortest path,100,Graphs.jl,7.291e-6,94 +min-cost flow,100,GraphsOptim.jl,0.00401075,7439.0 +min-cost flow,100,NetworkX,0.001786333043128252,7439.0 +shortest path,200,GraphsOptim.jl,0.049725583,80 +shortest path,200,Graphs.jl,1.6042e-5,80 +min-cost flow,200,GraphsOptim.jl,0.008109,8886.0 +min-cost flow,200,NetworkX,0.004978209035471082,8886.0 +shortest path,500,GraphsOptim.jl,0.122681625,104 +shortest path,500,Graphs.jl,4.6042e-5,104 +min-cost flow,500,GraphsOptim.jl,0.019684875,27758.0 +min-cost flow,500,NetworkX,0.024996375082992017,27758.0 +shortest path,1000,GraphsOptim.jl,0.2698078955,99 +shortest path,1000,Graphs.jl,0.000129417,99 +min-cost flow,1000,GraphsOptim.jl,0.0498085,45466.0 +min-cost flow,1000,NetworkX,0.05041183403227478,45466.0 +shortest path,2000,GraphsOptim.jl,0.5496791875,115 +shortest path,2000,Graphs.jl,0.000293375,115 +min-cost flow,2000,GraphsOptim.jl,0.131062833,99730.0 +min-cost flow,2000,NetworkX,0.1522331249434501,99730.0 +assignment,10,GraphsOptim.jl (JuMP),0.000908083,1.329931040873881 +assignment,10,GraphsOptim.jl (native),4.16e-7,1.329931040873881 +assignment,10,SciPy,1.166015863418579e-6,1.329931040873881 +assignment,20,GraphsOptim.jl (JuMP),0.00160425,1.3003183211649092 +assignment,20,GraphsOptim.jl (native),2.208e-6,1.3003183211649094 +assignment,20,SciPy,2.5830231606960297e-6,1.3003183211649094 +assignment,50,GraphsOptim.jl (JuMP),0.008272917,1.9305540673944215 +assignment,50,GraphsOptim.jl (native),1.2875e-5,1.9305540673944213 +assignment,50,SciPy,1.537497155368328e-5,1.9305540673944213 +assignment,100,GraphsOptim.jl (JuMP),0.03232875,1.5128939402479162 +assignment,100,GraphsOptim.jl (native),6.7375e-5,1.512893940247916 +assignment,100,SciPy,9.674998000264168e-5,1.512893940247916 +assignment,200,GraphsOptim.jl (JuMP),0.132599375,1.5860871489206083 +assignment,200,GraphsOptim.jl (native),0.000465416,1.5860871489206085 +assignment,200,SciPy,0.0005424169357866049,1.5860871489206085 +assignment,500,GraphsOptim.jl (JuMP),1.071073375,1.6176870888569141 +assignment,500,GraphsOptim.jl (native),0.004723708,1.617687088856912 +assignment,500,SciPy,0.004273458034731448,1.617687088856912 +assignment,1000,GraphsOptim.jl (JuMP),4.389844625,1.6825274033896138 +assignment,1000,GraphsOptim.jl (native),0.024060917,1.6825274033896118 +assignment,1000,SciPy,0.024344166042283177,1.6825274033896118 diff --git a/benchmark/run.jl b/benchmark/run.jl new file mode 100644 index 0000000..4de53dc --- /dev/null +++ b/benchmark/run.jl @@ -0,0 +1,47 @@ +# Benchmarks of GraphsOptim.jl against other Julia and Python libraries. +# +# Usage, from the root of the repository: +# +# julia --project=benchmark -e 'using Pkg; Pkg.instantiate()' +# julia --project=benchmark benchmark/run.jl [suite...] [--plots-only] +# +# Without arguments, all suites are run. The results are written to `benchmark/results/` +# and the figures to `docs/src/assets/benchmarks/`, both shown in the documentation. + +include("common.jl") +include("plots.jl") +include("suites/graph_matching.jl") +include("suites/edit_distance.jl") +include("suites/paths_flows.jl") +include("suites/cliques.jl") + +const SUITES = Dict( + "graph_matching" => (run_graph_matching, plot_graph_matching), + "edit_distance" => (run_edit_distance, plot_edit_distance), + "paths_flows" => (run_paths_flows, plot_paths_flows), + "cliques" => (run_cliques, plot_cliques), +) + +function main(args) + plots_only = "--plots-only" in args + names = filter(!=("--plots-only"), args) + isempty(names) && (names = sort(collect(keys(SUITES)))) + for name in names + haskey(SUITES, name) || error( + "Unknown suite $name, available suites: $(join(sort(collect(keys(SUITES))), ", "))", + ) + end + mkpath(FIGURES_DIR) + for name in names + run_suite, plot_suite = SUITES[name] + if !plots_only + @info "Running suite $name" + run_suite() + end + plot_suite() + end + plots_only || write_environment() + return nothing +end + +main(ARGS) diff --git a/benchmark/suites/cliques.jl b/benchmark/suites/cliques.jl new file mode 100644 index 0000000..30773ee --- /dev/null +++ b/benchmark/suites/cliques.jl @@ -0,0 +1,132 @@ +# Maximum weight clique, maximum independent set and minimum vertex cover: the integer +# programs of GraphsOptim.jl against exact algorithms and heuristics from Graphs.jl and NetworkX. + +is_clique(g, c) = all(has_edge(g, u, v) for u in c, v in c if u < v) +is_independent_set(g, s) = !any(has_edge(g, u, v) for u in s, v in s if u < v) +function is_vertex_cover(g, c) + return (cover = Set(c); all(src(e) in cover || dst(e) in cover for e in edges(g))) +end + +pyvertices(s) = pyconvert(Vector{Int}, pylist(s)) .+ 1 + +""" + clique_methods(g, w) + +Methods computing a maximum weight clique of `g` with vertex weights `w`, all exact. +Each returns a vector of vertices. +""" +function clique_methods(g, w) + nx = pyimport("networkx") + h = to_networkx(g; vertex_weights=w) + return [ + "GraphsOptim.jl" => + () -> maximum_weight_clique(g; vertex_weights=w, optimizer=OPTIMIZER), + "Graphs.jl (all maximal cliques)" => + () -> argmax(c -> sum(w[c]), maximal_cliques(g)), + "NetworkX" => () -> pyvertices(nx.max_weight_clique(h; weight="weight")[0]), + ] +end + +"Methods computing a (maximum or maximal) independent set, as pairs `name => (exact, f)`." +function independent_set_methods(g) + nx = pyimport("networkx") + h = to_networkx(g) + hc = nx.complement(h) + return [ + "GraphsOptim.jl" => + (true, () -> maximum_weight_independent_set(g; optimizer=OPTIMIZER)), + "NetworkX" => (true, () -> pyvertices(nx.max_weight_clique(hc; weight=nothing)[0])), + "Graphs.jl (greedy)" => (false, () -> independent_set(g, DegreeIndependentSet())), + "NetworkX (approximation)" => + (false, () -> pyvertices(nx.approximation.maximum_independent_set(h))), + ] +end + +"Methods computing a (minimum or small) vertex cover, as pairs `name => (exact, f)`." +function vertex_cover_methods(g) + nx = pyimport("networkx") + h = to_networkx(g) + return [ + "GraphsOptim.jl" => (true, () -> min_vertex_cover(g; optimizer=OPTIMIZER)), + "Graphs.jl (greedy)" => (false, () -> vertex_cover(g, DegreeVertexCover())), + "Graphs.jl (2-approximation)" => + (false, () -> vertex_cover(g, RandomVertexCover(); rng=Xoshiro(0))), + "NetworkX (2-approximation)" => + (false, () -> pyvertices(nx.approximation.min_weighted_vertex_cover(h))), + ] +end + +function run_methods!(rows, problem, n, methods, skipped; value, valid) + for (name, (exact, f)) in methods + name in skipped && continue + time, solution = measure(f) + too_slow(time) && push!(skipped, name) + ok = !isnothing(solution) + push!( + rows, + (; + problem, + n, + method=name, + exact, + time, + value=ok ? value(solution) : NaN, + valid=ok && valid(solution), + ), + ) + end +end + +function run_cliques() + # compilation + let g = erdos_renyi(10, 0.5; rng=Xoshiro(0)), w = ones(Int, 10) + foreach(((_, f),) -> f(), clique_methods(g, w)) + foreach(((_, (_, f)),) -> f(), independent_set_methods(g)) + foreach(((_, (_, f)),) -> f(), vertex_cover_methods(g)) + end + + rows = [] + skipped = Set{String}() + for n in [20, 40, 60, 80, 100, 150, 200] + rng = Xoshiro(n) + g = erdos_renyi(n, 0.5; rng) + w = rand(rng, 1:10, n) + methods = [name => (true, f) for (name, f) in clique_methods(g, w)] + run_methods!( + rows, + "maximum weight clique", + n, + methods, + skipped; + value=c -> sum(w[c]), + valid=c -> is_clique(g, c), + ) + @info "Clique" n + end + + skipped_is, skipped_vc = Set{String}(), Set{String}() + for n in [25, 50, 100, 200, 400, 800] + g = erdos_renyi(n, 5 / n; rng=Xoshiro(n)) + run_methods!( + rows, + "maximum independent set", + n, + independent_set_methods(g), + skipped_is; + value=length, + valid=s -> is_independent_set(g, s), + ) + run_methods!( + rows, + "minimum vertex cover", + n, + vertex_cover_methods(g), + skipped_vc; + value=length, + valid=c -> is_vertex_cover(g, c), + ) + @info "Independent set and vertex cover" n + end + write_csv("cliques.csv", rows) + return nothing +end diff --git a/benchmark/suites/edit_distance.jl b/benchmark/suites/edit_distance.jl new file mode 100644 index 0000000..6b480c8 --- /dev/null +++ b/benchmark/suites/edit_distance.jl @@ -0,0 +1,69 @@ +# Graph edit distance: the ILP formulations of GraphsOptim.jl against the exact +# algorithm of NetworkX, with uniform edit costs. + +const GED_FORMULATIONS = [F1, F1prime, F1plus, F2, F2minus, F2plus, FORI] + +""" + edited_graph(rng, g, k) + +Return a copy of `g` where `k` random vertex pairs have been flipped (edge added or removed), and the vertices shuffled. +""" +function edited_graph(rng::AbstractRNG, g::SimpleGraph, k::Integer) + h = copy(g) + pairs = [(i, j) for j in vertices(g) for i in 1:(j - 1)] + for (i, j) in pairs[randperm(rng, length(pairs))[1:k]] + has_edge(h, i, j) ? rem_edge!(h, i, j) : add_edge!(h, i, j) + end + return h[randperm(rng, nv(h))] +end + +function run_edit_distance() + nx = pyimport("networkx") + # compilation + let g = erdos_renyi(4, 0.5; rng=Xoshiro(0)) + for formulation in GED_FORMULATIONS + GraphsOptim.edit_distance( + g, edited_graph(Xoshiro(0), g, 1); formulation, optimizer=OPTIMIZER + ) + end + nx.graph_edit_distance(to_networkx(g), to_networkx(edited_graph(Xoshiro(0), g, 1))) + end + + rows = [] + timed_out = Set{String}() + for n in 4:2:24 + rng = Xoshiro(n) + G = erdos_renyi(n, 0.4; rng) + H = edited_graph(rng, G, max(2, n ÷ 3)) + for formulation in GED_FORMULATIONS + method = string(nameof(formulation)) + method in timed_out && continue + t, res = measure( + () -> GraphsOptim.edit_distance(G, H; formulation, optimizer=OPTIMIZER); + min_time=0, + ) + too_slow(t) && push!(timed_out, method) + value = isnothing(res) ? NaN : res.objective_value + push!(rows, (; n, edges_G=ne(G), edges_H=ne(H), method, time=t, value)) + end + if !("NetworkX" in timed_out) + t, value = pymeasure( + nx.graph_edit_distance, + to_networkx(G), + to_networkx(H); + timeout=TIME_LIMIT, + min_time=0, + ) + # On timeout, NetworkX returns the best value found so far, which is not proven optimal + value = t >= TIME_LIMIT ? NaN : pyconvert(Float64, value) + t = t >= TIME_LIMIT ? NaN : t + too_slow(t) && push!(timed_out, "NetworkX") + push!( + rows, (; n, edges_G=ne(G), edges_H=ne(H), method="NetworkX", time=t, value) + ) + end + @info "Edit distance" n + end + write_csv("edit_distance.csv", rows) + return nothing +end diff --git a/benchmark/suites/graph_matching.jl b/benchmark/suites/graph_matching.jl new file mode 100644 index 0000000..d880162 --- /dev/null +++ b/benchmark/suites/graph_matching.jl @@ -0,0 +1,80 @@ +# Graph matching: FAQ and GOAT against graspologic. + +""" + correlated_er(rng, n, p, ρ) + +Sample a pair of `ρ`-correlated Erdős–Rényi graphs with edge probability `p`, and shuffle the vertices of the second one. +Return the adjacency matrices `A`, `B` and the ground truth `σ` such that vertex `i` of `A` corresponds to vertex `σ[i]` of `B`. +""" +function correlated_er(rng::AbstractRNG, n::Integer, p::Real, ρ::Real) + A, B = zeros(n, n), zeros(n, n) + for j in 1:n, i in 1:(j - 1) + a = rand(rng) < p + b = a ? rand(rng) < p + ρ * (1 - p) : rand(rng) < p * (1 - ρ) + A[i, j] = A[j, i] = a + B[i, j] = B[j, i] = b + end + π = randperm(rng, n) + return A, B[π, π], invperm(π) +end + +"Number of edges present in exactly one of the two aligned graphs, i.e. `‖A - P B Pᵀ‖² / 2`." +disagreements(A, B, σ) = round(Int, sum(abs2, A - B[σ, σ]) / 2) + +function graph_matching_rows(algo, A, B, truth; seed, min_time=1.0) + name = string(nameof(typeof(algo))) + graspologic = pyimport("graspologic.match") + np = pyimport("numpy") + t_jl, (P, _, _) = measure(() -> graph_matching(algo, A, B); min_time) + σ_jl = [argmax(view(P, i, :)) for i in axes(P, 1)] + t_py, res = pymeasure( + graspologic.graph_match, + np.ascontiguousarray(A), + np.ascontiguousarray(B); + transport=algo isa GOAT, + rng=seed, + min_time, + ) + σ_py = pyconvert(Vector{Int}, res.indices_B[np.argsort(res.indices_A)]) .+ 1 + return [ + (; + algorithm=name, + library, + time, + accuracy=mean(σ .== truth), + disagreements=disagreements(A, B, σ), + disagreements_truth=disagreements(A, B, truth), + ) for + (library, time, σ) in (("GraphsOptim.jl", t_jl, σ_jl), ("graspologic", t_py, σ_py)) + ] +end + +function run_graph_matching() + # compilation + let (A, B, truth) = correlated_er(Xoshiro(0), 20, 0.3, 0.9) + graph_matching_rows(FAQ(), A, B, truth; seed=0, min_time=0) + graph_matching_rows(GOAT(), A, B, truth; seed=0, min_time=0) + end + + # Running time as the number of vertices grows + scaling = [] + for n in [50, 100, 200, 500, 1000], algo in (FAQ(), GOAT()) + A, B, truth = correlated_er(Xoshiro(n), n, 0.1, 0.9) + for row in graph_matching_rows(algo, A, B, truth; seed=n) + push!(scaling, (; n, row...)) + end + @info "Graph matching" n algo + end + write_csv("graph_matching_scaling.csv", scaling) + + # Matching quality as the correlation between the graphs varies + quality = [] + for ρ in 0.5:0.1:1.0, seed in 1:10, algo in (FAQ(), GOAT()) + A, B, truth = correlated_er(Xoshiro(seed), 150, 0.1, ρ) + for row in graph_matching_rows(algo, A, B, truth; seed, min_time=0) + push!(quality, (; ρ, seed, row...)) + end + end + write_csv("graph_matching_quality.csv", quality) + return nothing +end diff --git a/benchmark/suites/paths_flows.jl b/benchmark/suites/paths_flows.jl new file mode 100644 index 0000000..6c96656 --- /dev/null +++ b/benchmark/suites/paths_flows.jl @@ -0,0 +1,151 @@ +# Shortest path, min-cost flow and assignment: the linear programs of GraphsOptim.jl +# against combinatorial algorithms from Graphs.jl, NetworkX and SciPy. + +"Random directed graph with average out-degree `d`, plus a Hamiltonian cycle so that it is strongly connected." +function random_digraph(rng::AbstractRNG, n::Integer, d::Real) + g = erdos_renyi(n, d / n; is_directed=true, rng) + for i in 1:n + add_edge!(g, i, mod1(i + 1, n)) + end + return g +end + +"Random integer edge weights in `1:wmax`, as a dense matrix (zero outside of the edges)." +function random_edge_weights(rng::AbstractRNG, g::AbstractGraph, wmax::Integer) + W = zeros(Int, nv(g), nv(g)) + for e in edges(g) + W[src(e), dst(e)] = rand(rng, 1:wmax) + end + return W +end + +path_cost(path, W) = sum(W[path[k], path[k + 1]] for k in 1:(length(path) - 1)) + +function shortest_path_rows(n) + rng = Xoshiro(n) + g = random_digraph(rng, n, 5) + W = random_edge_weights(rng, g, 100) + s, t = 1, n ÷ 2 + rows = [] + time, path = measure(() -> shortest_path(g, s, t, W; optimizer=OPTIMIZER)) + push!( + rows, + (; + problem="shortest path", + n, + library="GraphsOptim.jl", + time, + value=path_cost(path, W), + ), + ) + time, path = measure(() -> enumerate_paths(dijkstra_shortest_paths(g, s, W), t)) + push!( + rows, + (; problem="shortest path", n, library="Graphs.jl", time, value=path_cost(path, W)), + ) + return rows +end + +function min_cost_flow_rows(n) + rng = Xoshiro(n) + g = random_digraph(rng, n, 5) + cost = random_edge_weights(rng, g, 100) + capacity = random_edge_weights(rng, g, 20) + # a tenth of the vertices are sources, another tenth are sinks + demand = zeros(Int, n) + vs = randperm(rng, n) + k = max(1, n ÷ 10) + for i in 1:k + supply = rand(rng, 1:10) + demand[vs[i]] -= supply + demand[vs[k + i]] += supply + end + # the Hamiltonian cycle can carry all the flow, at a high cost, so the problem is feasible + for i in 1:n + capacity[i, mod1(i + 1, n)] = sum(abs, demand) + cost[i, mod1(i + 1, n)] = 100 + end + rows = [] + time, flow = measure( + () -> + min_cost_flow(g, demand, cost, zeros(Int, n, n), capacity; optimizer=OPTIMIZER), + ) + value = isnothing(flow) ? NaN : sum(cost .* flow) + push!(rows, (; problem="min-cost flow", n, library="GraphsOptim.jl", time, value)) + + nx = pyimport("networkx") + h = nx.DiGraph() + for v in 1:n + h.add_node(v - 1; demand=demand[v]) + end + for e in edges(g) + i, j = src(e), dst(e) + h.add_edge(i - 1, j - 1; weight=cost[i, j], capacity=capacity[i, j]) + end + time, res = pymeasure(nx.network_simplex, h) + push!( + rows, + (; + problem="min-cost flow", + n, + library="NetworkX", + time, + value=pyconvert(Float64, res[0]), + ), + ) + return rows +end + +function assignment_rows(n; with_jump) + C = rand(Xoshiro(n), n, n) + rows = [] + if with_jump + time, P = measure(() -> min_cost_assignment(C; integer=false, optimizer=OPTIMIZER)) + value = isnothing(P) ? NaN : dot(C, P) + push!( + rows, (; problem="assignment", n, library="GraphsOptim.jl (JuMP)", time, value) + ) + end + time, σ = measure(() -> GraphsOptim.linear_sum_assignment(C)) + value = sum(C[i, σ[i]] for i in 1:n) + push!(rows, (; problem="assignment", n, library="GraphsOptim.jl (native)", time, value)) + scipy = pyimport("scipy.optimize") + np = pyimport("numpy") + time, res = pymeasure(scipy.linear_sum_assignment, np.ascontiguousarray(C)) + σ = pyconvert(Vector{Int}, res[1]) .+ 1 + push!( + rows, + (; + problem="assignment", + n, + library="SciPy", + time, + value=sum(C[i, σ[i]] for i in 1:n), + ), + ) + return rows +end + +function run_paths_flows() + # compilation + shortest_path_rows(10) + min_cost_flow_rows(10) + assignment_rows(10; with_jump=true) + + rows = [] + for n in [100, 200, 500, 1000, 2000] + append!(rows, shortest_path_rows(n)) + append!(rows, min_cost_flow_rows(n)) + @info "Paths and flows" n + end + jump_ok = true + for n in [10, 20, 50, 100, 200, 500, 1000] + new_rows = assignment_rows(n; with_jump=jump_ok) + # stop using JuMP once it hits the time limit + jump_ok = jump_ok && !isnan(first(new_rows).time) + append!(rows, new_rows) + @info "Assignment" n + end + write_csv("paths_flows.csv", rows) + return nothing +end diff --git a/docs/make.jl b/docs/make.jl index d932cc2..f7a504d 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -31,7 +31,11 @@ makedocs(; prettyurls=get(ENV, "CI", "false") == "true", canonical="https://juliagraphs.org/GraphsOptim.jl", ), - pages=["Home" => "index.md", "Algorithms" => "algorithms.md"], + pages=[ + "Home" => "index.md", + "Algorithms" => "algorithms.md", + "Benchmarks" => "benchmarks.md", + ], ) deploydocs(; diff --git a/docs/src/algorithms.md b/docs/src/algorithms.md index 0daef5f..b43729e 100644 --- a/docs/src/algorithms.md +++ b/docs/src/algorithms.md @@ -59,6 +59,7 @@ A special case of minimum cost flow without edge capacities, and where vertex de ```@docs min_cost_assignment GraphsOptim.min_cost_assignment! +GraphsOptim.linear_sum_assignment ``` ## Minimum Vertex Cover diff --git a/docs/src/assets/benchmarks/cliques_quality.svg b/docs/src/assets/benchmarks/cliques_quality.svg new file mode 100644 index 0000000..062fc61 --- /dev/null +++ b/docs/src/assets/benchmarks/cliques_quality.svg @@ -0,0 +1,898 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/assets/benchmarks/cliques_time.svg b/docs/src/assets/benchmarks/cliques_time.svg new file mode 100644 index 0000000..dab0201 --- /dev/null +++ b/docs/src/assets/benchmarks/cliques_time.svg @@ -0,0 +1,1100 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/assets/benchmarks/edit_distance_time.svg b/docs/src/assets/benchmarks/edit_distance_time.svg new file mode 100644 index 0000000..e7c429b --- /dev/null +++ b/docs/src/assets/benchmarks/edit_distance_time.svg @@ -0,0 +1,390 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/assets/benchmarks/graph_matching_accuracy.svg b/docs/src/assets/benchmarks/graph_matching_accuracy.svg new file mode 100644 index 0000000..2e3dcf5 --- /dev/null +++ b/docs/src/assets/benchmarks/graph_matching_accuracy.svg @@ -0,0 +1,406 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/assets/benchmarks/graph_matching_time.svg b/docs/src/assets/benchmarks/graph_matching_time.svg new file mode 100644 index 0000000..fe7043b --- /dev/null +++ b/docs/src/assets/benchmarks/graph_matching_time.svg @@ -0,0 +1,447 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/assets/benchmarks/paths_flows_time.svg b/docs/src/assets/benchmarks/paths_flows_time.svg new file mode 100644 index 0000000..68d1930 --- /dev/null +++ b/docs/src/assets/benchmarks/paths_flows_time.svg @@ -0,0 +1,848 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/benchmarks.md b/docs/src/benchmarks.md new file mode 100644 index 0000000..3e13430 --- /dev/null +++ b/docs/src/benchmarks.md @@ -0,0 +1,138 @@ +# Benchmarks + +GraphsOptim.jl solves graph problems by formulating them as mathematical programs. +These benchmarks compare its algorithms with dedicated implementations from other Julia and Python libraries, in terms of running time and, when the reference is a heuristic, of solution quality. + +## Where to find them + +The benchmarks live in the [`benchmark/`](https://github.com/JuliaGraphs/GraphsOptim.jl/tree/main/benchmark) folder of the repository, which has its own Julia environment: + +| File | Content | +| :--- | :--- | +| `benchmark/run.jl` | entry point, runs the suites and draws the figures | +| `benchmark/suites/*.jl` | one suite per family of problems: instances, methods and metrics | +| `benchmark/common.jl` | shared tools: timing, calls to Python, CSV files | +| `benchmark/plots.jl` | figures, saved in `docs/src/assets/benchmarks/` | +| `benchmark/results/*.csv` | raw results displayed on this page | + +To run all suites, or only some of them, from the root of the repository: + +```bash +julia --project=benchmark -e 'using Pkg; Pkg.instantiate()' +julia --project=benchmark benchmark/run.jl # all suites +julia --project=benchmark benchmark/run.jl cliques # a single suite +julia --project=benchmark benchmark/run.jl --plots-only # only redraw the figures +``` + +The Python libraries ([graspologic](https://graspologic-org.github.io/graspologic/), [NetworkX](https://networkx.org/) and [SciPy](https://scipy.org/)) are installed automatically by [CondaPkg.jl](https://github.com/JuliaPy/CondaPkg.jl) and called through [PythonCall.jl](https://github.com/JuliaPy/PythonCall.jl). + +### Methodology + +- All methods are run on the same random instances, with their default parameters. +- Running times are medians over repeated runs (at least one second or seven runs), measured after compilation. Python functions are timed on the Python side, so that the conversion of their inputs is not counted. +- All libraries use a single BLAS thread. +- The optimizer of GraphsOptim.jl is [HiGHS](https://highs.dev/), with a time limit of 60 seconds. A method that fails or takes more than 6 seconds on an instance is not run on the larger ones, hence the missing points in the figures. + +The current results were obtained with: + +```@eval +using GraphsOptim, Markdown +Markdown.parse(read(joinpath(pkgdir(GraphsOptim), "benchmark", "results", "environment.md"), String)) +``` + +```@setup bench +using GraphsOptim, Markdown + +function read_results(file) + lines = readlines(joinpath(pkgdir(GraphsOptim), "benchmark", "results", file)) + header = Tuple(Symbol.(split(lines[1], ","))) + parsevalue(x) = something(tryparse(Int, x), tryparse(Float64, x), String(x)) + return [NamedTuple{header}(Tuple(parsevalue.(split(line, ",")))) for line in lines[2:end]] +end + +""" +Check that the methods in `rows` found the same value on each instance (grouped by `instance`). +""" +function agreement(rows; instance, what, who="all methods") + mismatches = String[] + for key in unique(instance.(rows)) + values = [r.value for r in rows if instance(r) == key && !(r.value isa Float64 && isnan(r.value))] + if !all(v -> isapprox(v, first(values); rtol=1e-6), values) + push!(mismatches, string(key)) + end + end + if isempty(mismatches) + return Markdown.parse("On every instance, $who found the same $what.") + else + return Markdown.parse("!!! warning\n The methods disagree on the $what for: $(join(mismatches, ", ")).") + end +end +``` + +## Graph matching + +[`graph_matching`](@ref) with `FAQ()` and `GOAT()`, against `graspologic.match.graph_match` (with `transport=True` for GOAT). +The instances are pairs of correlated Erdős–Rényi graphs with edge probability ``p = 0.1``: the second graph keeps each edge of the first one with a probability depending on the correlation ``\rho``, and its vertices are shuffled. + +Running time for ``\rho = 0.9``: + +![Running time of graph matching](assets/benchmarks/graph_matching_time.svg) + +Fraction of vertices matched to their ground-truth counterpart, for ``n = 150`` vertices (mean over 10 instances): + +![Accuracy of graph matching](assets/benchmarks/graph_matching_accuracy.svg) + +FAQ, initialized at the flat doubly stochastic matrix, only recovers the alignment of identical graphs; GOAT recovers it as soon as the correlation is high enough. Both libraries return alignments of the same quality. + +## Graph edit distance + +[`GraphsOptim.edit_distance`](@ref) with each of its integer programming formulations, against the exact algorithm `networkx.graph_edit_distance`, with uniform edit costs. +The instances are Erdős–Rényi graphs with edge probability ``0.4``, and copies of them where ``n / 3`` random vertex pairs have been flipped (an edge is added or removed). + +![Running time of graph edit distance](assets/benchmarks/edit_distance_time.svg) + +```@example bench +agreement(read_results("edit_distance.csv"); instance=r -> r.n, what="edit distance") +``` + +## Shortest path, min-cost flow and assignment + +- [`shortest_path`](@ref) against Dijkstra's algorithm from Graphs.jl (`dijkstra_shortest_paths`), on random directed graphs with average degree 5 and integer edge costs. +- [`min_cost_flow`](@ref) against the network simplex of NetworkX (`networkx.network_simplex`), on the same kind of graphs with integer capacities, where a tenth of the vertices are sources and another tenth are sinks. +- [`min_cost_assignment`](@ref), which solves a linear program with JuMP, and [`GraphsOptim.linear_sum_assignment`](@ref), against `scipy.optimize.linear_sum_assignment`, on uniformly random cost matrices. + +![Running time of shortest path, min-cost flow and assignment](assets/benchmarks/paths_flows_time.svg) + +```@example bench +agreement(read_results("paths_flows.csv"); instance=r -> (r.problem, r.n), what="optimal cost") +``` + +## Cliques, independent sets and vertex covers + +- [`maximum_weight_clique`](@ref) against the exact algorithms `networkx.max_weight_clique` and the enumeration of all maximal cliques with Graphs.jl (`maximal_cliques`), on Erdős–Rényi graphs with edge probability ``0.5`` and integer vertex weights. +- [`maximum_weight_independent_set`](@ref) (with unit weights) and [`min_vertex_cover`](@ref) against exact and approximate algorithms from NetworkX and the heuristics of Graphs.jl, on Erdős–Rényi graphs with average degree 5. + +![Running time of clique, independent set and vertex cover](assets/benchmarks/cliques_time.svg) + +```@example bench +cliques = read_results("cliques.csv") +agreement( + filter(r -> r.exact == "true", cliques); + instance=r -> (r.problem, r.n), + what="optimal value", + who="all exact methods", +) +``` + +```@example bench +invalid = [(r.problem, r.n, r.method) for r in cliques if r.valid != "true" && !(r.value isa Float64 && isnan(r.value))] +if isempty(invalid) + Markdown.parse("All returned solutions are valid cliques, independent sets and vertex covers.") +else + Markdown.parse("!!! warning\n Invalid solutions: $(join(string.(invalid), ", ")).") +end +``` + +The integer programs of GraphsOptim.jl are exact, but the heuristics are much faster. Their solutions are compared to the optimum below (higher is better for independent sets, lower is better for vertex covers): + +![Quality of the independent set and vertex cover heuristics](assets/benchmarks/cliques_quality.svg) diff --git a/src/GraphsOptim.jl b/src/GraphsOptim.jl index 02496fd..b645ca6 100644 --- a/src/GraphsOptim.jl +++ b/src/GraphsOptim.jl @@ -17,10 +17,11 @@ using JuMP: set_silent, optimize!, termination_status, value using JuMP: set_optimizer, objective_value using JuMP: @variable, @constraint, @objective using JuMP: VariableRef -using JuMP: num_variables -using LinearAlgebra: norm, tr, dot +using JuMP: num_variables, backend, index +using LinearAlgebra: norm, dot, mul!, issymmetric using MathOptInterface: OPTIMAL -using SparseArrays: sparse +import MathOptInterface as MOI +using SparseArrays: sparse, SparseMatrixCSC, rowvals, nonzeros, nzrange using OptimalTransport: sinkhorn export min_cost_flow diff --git a/src/assignment.jl b/src/assignment.jl index 9e50d53..e5dbf94 100644 --- a/src/assignment.jl +++ b/src/assignment.jl @@ -13,10 +13,15 @@ function min_cost_assignment!( model::Model, edge_cost::AbstractMatrix; var_name, integer::Bool ) n, m = size(edge_cost) - x = @variable(model, [i in 1:n, j in 1:m]; integer=integer, base_name=String(var_name)) + x = @variable( + model, + [i in 1:n, j in 1:m]; + lower_bound=0, + integer=integer, + base_name=String(var_name) + ) model[Symbol(var_name)] = x @objective(model, Min, dot(edge_cost, x)) - @constraint(model, x .>= 0) @constraint(model, [i in 1:n], sum(x[i, :]) == 1) @constraint(model, [j in 1:m], sum(x[:, j]) == 1) return nothing @@ -49,3 +54,95 @@ function min_cost_assignment( @assert termination_status(model) == OPTIMAL return value.(model[:assignment]) end + +""" + linear_sum_assignment(cost; maximize=false) + +Solve the square linear assignment problem with a shortest augmenting path algorithm, without relying on a JuMP optimizer. + +Return a vector `σ` such that row `i` is assigned to column `σ[i]`, minimizing (or maximizing) `∑ᵢ cost[i, σ[i]]`. + +# References + +- Crouse, "On implementing 2D rectangular assignment algorithms", IEEE Transactions on Aerospace and Electronic Systems (2016), also used by `scipy.optimize.linear_sum_assignment` +""" +function linear_sum_assignment(cost::AbstractMatrix{<:Real}; maximize::Bool=false) + n = size(cost, 1) + size(cost, 2) == n || throw(DimensionMismatch("the cost matrix must be square")) + T = float(eltype(cost)) + s = maximize ? -one(T) : one(T) + # The columns of `cost` play the role of workers and its rows the role of jobs, + # so that the innermost loop runs down a column (memory-contiguous in Julia). + u = zeros(T, n) # dual variables of the workers + v = zeros(T, n) # dual variables of the jobs + shortest = Vector{T}(undef, n) + path = zeros(Int, n) + job4worker = zeros(Int, n) + worker4job = zeros(Int, n) + visited_workers = Vector{Bool}(undef, n) + visited_jobs = Vector{Bool}(undef, n) + remaining = Vector{Int}(undef, n) + @inbounds for current in 1:n + fill!(shortest, typemax(T)) + fill!(visited_workers, false) + fill!(visited_jobs, false) + for k in 1:n + remaining[k] = n - k + 1 + end + num_remaining = n + min_val = zero(T) + w = current + sink = 0 + # Find a shortest augmenting path starting from worker `current` + while sink == 0 + visited_workers[w] = true + index = 0 + lowest = typemax(T) + uw = u[w] + for k in 1:num_remaining + j = remaining[k] + r = min_val + s * cost[j, w] - uw - v[j] + if r < shortest[j] + path[j] = w + shortest[j] = r + end + if shortest[j] < lowest || (shortest[j] == lowest && worker4job[j] == 0) + lowest = shortest[j] + index = k + end + end + min_val = lowest + isfinite(min_val) || throw(ArgumentError("the cost matrix is infeasible")) + j = remaining[index] + if worker4job[j] == 0 + sink = j + else + w = worker4job[j] + end + visited_jobs[j] = true + remaining[index] = remaining[num_remaining] + num_remaining -= 1 + end + # Update the dual variables + u[current] += min_val + for w in 1:n + if visited_workers[w] && w != current + u[w] += min_val - shortest[job4worker[w]] + end + end + for j in 1:n + if visited_jobs[j] + v[j] -= min_val - shortest[j] + end + end + # Augment the assignment along the path + j = sink + while true + w = path[j] + worker4job[j] = w + job4worker[w], j = j, job4worker[w] + w == current && break + end + end + return worker4job +end diff --git a/src/graph_matching.jl b/src/graph_matching.jl index 75163f3..a04b69d 100644 --- a/src/graph_matching.jl +++ b/src/graph_matching.jl @@ -4,56 +4,266 @@ struct FAQ <: GraphMatchingAlgorithm end struct GOAT <: GraphMatchingAlgorithm end +""" + graph_matching_step_size(fP, fQ, cPQ) + +Return the step size `α ∈ [0, 1]` maximizing the quadratic objective `f(αP + (1-α)Q) = α² fP + (1-α)² fQ + α(1-α) cPQ` of graph matching along the Frank-Wolfe segment, where `fP = f(P)`, `fQ = f(Q)` and `cPQ = ⟨∇f(P), Q⟩`. +""" +function graph_matching_step_size(fP::Real, fQ::Real, cPQ::Real) + a = fP + fQ - cPQ + b = cPQ - 2fQ + if a < 0 + # concave parabola: clamp its vertex to the segment + return clamp(-b / (2a), zero(a), one(a)) + else + # convex or linear: the best endpoint wins + return fP > fQ ? one(a) : zero(a) + end +end + """ graph_matching_step_size(A, B, P, Q) -Given the adjacency matrices `A` and `B`, the doubly stochastic matrix `P` and the direction matrix `Q`, return the step size of the Frank-Wolfe method for graph matching algorithms. +Given the adjacency matrices `A` and `B`, the doubly stochastic matrix `P` and the direction matrix `Q`, return the step size of the Frank-Wolfe method for graph matching algorithms, i.e. the `α ∈ [0, 1]` maximizing `tr(A X Bᵀ Xᵀ)` for `X = αP + (1-α)Q`. """ function graph_matching_step_size( A::AbstractMatrix, B::AbstractMatrix, P::AbstractMatrix, Q::AbstractMatrix ) - R = P - Q - a = -tr(A * R * B' * R') - b = -tr(A * Q * B' * R' + A * R * B' * Q') - if a > eps(typeof(a)) - α = min(max(b / (2a), 0.0), 1.0) - elseif a + b > eps(typeof(a + b)) - α = 0.0 + ∇fP = A * P * B' + A' * P * B + fP = dot(∇fP, P) / 2 + fQ = dot(A * Q * B', Q) + cPQ = dot(∇fP, Q) + return graph_matching_step_size(fP, fQ, cPQ) +end + +_dense_or_sparse(M::SparseMatrixCSC, ::Type{T}) where {T} = SparseMatrixCSC{T}(M) +_dense_or_sparse(M::AbstractMatrix, ::Type{T}) where {T} = convert(Matrix{T}, M) + +""" + _permute_rows!(X, M, σ) + +Fill the dense matrix `X` with `M[σ, :]`. +""" +function _permute_rows!(X::AbstractMatrix, M::AbstractMatrix, σ::AbstractVector{<:Integer}) + @inbounds for j in axes(M, 2), i in axes(M, 1) + X[i, j] = M[σ[i], j] + end + return X +end + +function _permute_rows!(X::AbstractMatrix, M::SparseMatrixCSC, σ::AbstractVector{<:Integer}) + fill!(X, zero(eltype(X))) + σinv = invperm(σ) + rows, vals = rowvals(M), nonzeros(M) + @inbounds for j in axes(M, 2), k in nzrange(M, j) + X[σinv[rows[k]], j] = vals[k] + end + return X +end + +""" + _gradient!(G, X, A, B, Bt, P, symmetric) + +Overwrite `G` with the gradient `A P Bᵀ + Aᵀ P B` of `P ↦ tr(A P Bᵀ Pᵀ)`, using `X` as a buffer. +""" +function _gradient!(G, X, A, B, Bt, P::AbstractMatrix, symmetric::Bool) + mul!(X, P, Bt) + mul!(G, A, X) + if symmetric + G .*= 2 + else + mul!(X, P, B) + mul!(G, transpose(A), X, true, true) + end + return G +end + +# For a permutation `Q` with `Q[i, σ[i]] = 1`, we have `Q Bᵀ = Bᵀ[σ, :]` and `Q B = B[σ, :]` +function _gradient!(G, X, A, B, Bt, σ::AbstractVector{<:Integer}, symmetric::Bool) + _permute_rows!(X, Bt, σ) + mul!(G, A, X) + if symmetric + G .*= 2 else - α = 1.0 + _permute_rows!(X, B, σ) + mul!(G, transpose(A), X, true, true) end - return α + return G end -function graph_matching_search_direction(::FAQ, ∇f; optimizer, kwargs...) - return min_cost_assignment(∇f; integer=false, optimizer=optimizer) +# Barycenter `J = 𝟏𝟏ᵀ / n`: the gradient has rank two and needs no matrix product +function _flat_gradient!(G, A, B) + n = size(A, 1) + rA, rB = Vector(vec(sum(A; dims=2))), Vector(vec(sum(B; dims=2))) + cA, cB = Vector(vec(sum(A; dims=1))), Vector(vec(sum(B; dims=1))) + G .= (rA .* rB' .+ cA .* cB') ./ n + return G end -function graph_matching_search_direction( - ::GOAT, ∇f::AbstractMatrix{T}; regularization, max_iter_sinkhorn, kwargs... -) where {T} - m = size(∇f, 1) - return sinkhorn( - ones(T, m), - ones(T, m), - -∇f, - -inv(regularization / maximum(abs, ∇f)); - maxiter=max_iter_sinkhorn, - ) +""" + _assignment_model(n, optimizer) + +Build once the JuMP model of the linear assignment problem (a linear program over the doubly stochastic matrices), whose objective is then updated at every Frank-Wolfe iteration. +""" +function _assignment_model(n::Integer, optimizer) + model = Model(optimizer) + set_silent(model) + min_cost_assignment!(model, zeros(n, n); var_name=:assignment, integer=false) + return model +end + +""" + _argmax_permutation(score) + +Return the row-wise maxima `σ` of `score` if they form a permutation, and `nothing` otherwise. +In the first case, `σ` solves the assignment problem `max ∑ᵢ score[i, σ[i]]`, since the sum of the row-wise maxima is an upper bound. +""" +function _argmax_permutation(score::AbstractMatrix) + σ = [argmax(view(score, i, :)) for i in axes(score, 1)] + return isperm(σ) ? σ : nothing +end + +""" + _max_assignment!(model, score) + +Return an optimal solution of `max ⟨score, X⟩` over the doubly stochastic matrices `X`. + +If `model` is `nothing`, the problem is solved with [`GraphsOptim.linear_sum_assignment`](@ref). +Otherwise, the objective of the JuMP assignment `model` is updated before solving it again, which allows the optimizer to warm start from the previous solution. +""" +function _max_assignment!(model::Union{Nothing,Model}, score::AbstractMatrix) + σ = _argmax_permutation(score) + if isnothing(σ) && isnothing(model) + σ = linear_sum_assignment(score; maximize=true) + end + if !isnothing(σ) + Q = zeros(eltype(score), size(score)) + for i in axes(score, 1) + Q[i, σ[i]] = 1 + end + return Q + end + x = vec(model[:assignment]) + # Replacing the whole objective at once is much faster than modifying each coefficient + f = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(-vec(score), index.(x)), 0.0) + MOI.set(backend(model), MOI.ObjectiveFunction{typeof(f)}(), f) + optimize!(model) + @assert termination_status(model) == OPTIMAL + return reshape(value.(x), size(score)) +end + +""" + _permutation_vector(Q) + +Return `σ` such that `Q[i, σ[i]] = 1` if `Q` is (numerically) a permutation matrix, and `nothing` otherwise. +""" +function _permutation_vector(Q::AbstractMatrix) + σ = [argmax(view(Q, i, :)) for i in axes(Q, 1)] + isperm(σ) || return nothing + for i in axes(Q, 1) + isapprox(Q[i, σ[i]], 1; atol=1e-6) || return nothing + end + return σ +end + +""" + _distance(P, Q) + +Frobenius distance `‖P - Q‖` between two matrices, without allocating `P - Q`. +""" +function _distance(P::AbstractMatrix, Q::AbstractMatrix) + d = zero(promote_type(eltype(P), eltype(Q))) + for k in eachindex(P, Q) + d += abs2(P[k] - Q[k]) + end + return sqrt(d) +end + +function _step_cache(::FAQ, n, optimizer) + return isnothing(optimizer) ? nothing : _assignment_model(n, optimizer) +end +_step_cache(::GOAT, n, optimizer) = nothing + +""" + _frank_wolfe_step!(algo, P, G, GQ, X, A, B, Bt, symmetric, cache; kwargs...) + +Perform one Frank-Wolfe step: compute the search direction `Q` from the gradient `G = ∇f(P)`, store `∇f(Q)` in `GQ`, the optimal step size `α`, and overwrite `P` with `αP + (1-α)Q`. +Return `α` and `‖P - Q‖`. +""" +function _frank_wolfe_step!(::FAQ, P, G, GQ, X, A, B, Bt, symmetric, model; kwargs...) + n = size(P, 1) + # Search direction: a vertex Q of the Birkhoff polytope maximizing ⟨∇f(P), Q⟩ + Q = _max_assignment!(model, G) + σ = _permutation_vector(Q) + if isnothing(σ) + _gradient!(GQ, X, A, B, Bt, Q, symmetric) + else + # cheaper gradient: Q B and Q Bᵀ are row permutations + _gradient!(GQ, X, A, B, Bt, σ, symmetric) + end + fP = dot(G, P) / 2 + fQ = dot(GQ, Q) / 2 + cPQ = dot(G, Q) + dist_PQ = _distance(P, Q) + α = graph_matching_step_size(fP, fQ, cPQ) + P .= α .* P .+ (1 - α) .* Q + return α, dist_PQ +end + +function _frank_wolfe_step!( + ::GOAT, + P, + G, + GQ, + X, + A, + B, + Bt, + symmetric, + cache; + regularization, + max_iter_sinkhorn, + tol_sinkhorn, + kwargs..., +) + n = size(P, 1) + # Search direction: the entropic transport plan maximizing ⟨∇f(P), Q⟩, + # computed with the cost `max(G) - G ≥ 0` (shifting the cost does not change the plan) + m = maximum(abs, G) + if m > 0 + Q = sinkhorn( + ones(eltype(P), n), + ones(eltype(P), n), + maximum(G) .- G, + m / regularization; + atol=tol_sinkhorn, + maxiter=max_iter_sinkhorn, + ) + else + Q = fill(one(eltype(P)) / n, n, n) + end + _gradient!(GQ, X, A, B, Bt, Q, symmetric) + fP = dot(G, P) / 2 + fQ = dot(GQ, Q) / 2 + cPQ = dot(G, Q) + dist_PQ = _distance(P, Q) + α = graph_matching_step_size(fP, fQ, cPQ) + P .= α .* P .+ (1 - α) .* Q + return α, dist_PQ end """ graph_matching( algo, A, B; optimizer, P_init, max_iter, tol, - max_iter_sinkhorn, regularizer + max_iter_sinkhorn, tol_sinkhorn, regularization ) -Compute an approximately optimal alignment between two graphs using one of two variants of Frank-Wolfe (FAQ or GOAT) +Compute an approximately optimal alignment between two graphs using one of two variants of Frank-Wolfe (FAQ or GOAT), by maximizing `tr(A P Bᵀ Pᵀ)` over the doubly stochastic matrices `P`. The output is a tuple that contains: - the permutation matrix `P` defining the alignment; -- the distance between the permuted graphs; +- the distance `‖A P - P B‖` between the permuted graphs; - a boolean indicating if the algorithm converged. # Arguments @@ -66,15 +276,16 @@ The output is a tuple that contains: For both algorithms: -- `optimizer`: JuMP-compatible solver (default is `HiGHS.Optimizer`) -- `P_init`: initialization matrix (default value is a flat doubly stochastic matrix). +- `optimizer`: JuMP-compatible solver for the linear assignment problems (default is `nothing`, which uses the built-in [`GraphsOptim.linear_sum_assignment`](@ref), much faster for this problem). +- `P_init`: initialization matrix (default value is the flat doubly stochastic matrix). - `max_iter`: maximum iterations of the Frank-Wolfe method (default value is 30). -- `tol`: tolerance for the convergence (default value is 0.1). +- `tol`: the algorithm stops when `‖Pₖ₊₁ - Pₖ‖ / √n < tol` (default value is 0.01). Only for GOAT: - `regularization`: penalty coefficient in the Sinkhorn algorithm (default value is 100.0). -- `max_iter_sinkhorn`: maximum iterations of the Sinkhorn algorithm (default value is 500). +- `max_iter_sinkhorn`: maximum iterations of the Sinkhorn algorithm (default value is 1000). +- `tol_sinkhorn`: tolerance on the marginals in the Sinkhorn algorithm in the 1-norm (default value is 0.05). # References @@ -85,36 +296,78 @@ function graph_matching( algo::GraphMatchingAlgorithm, A::AbstractMatrix, B::AbstractMatrix; - optimizer=HiGHS.Optimizer, - P_init::AbstractMatrix=flat_doubly_stochastic(size(A, 1)), + optimizer=nothing, + P_init::Union{Nothing,AbstractMatrix}=nothing, max_iter::Integer=30, - tol::Real=0.1, + tol::Real=0.01, regularization::Real=100.0, - max_iter_sinkhorn::Integer=500, + max_iter_sinkhorn::Integer=1000, + tol_sinkhorn::Real=0.05, ) - converged = false + n = size(A, 1) + size(A) == size(B) == (n, n) || + throw(DimensionMismatch("A and B must be square matrices of the same size")) + T = float(promote_type(eltype(A), eltype(B))) + A = _dense_or_sparse(A, T) + B = _dense_or_sparse(B, T) + symmetric = issymmetric(A) && issymmetric(B) + Bt = symmetric ? B : copy(transpose(B)) + G = Matrix{T}(undef, n, n) # gradient at P + GQ = Matrix{T}(undef, n, n) # gradient at Q + X = Matrix{T}(undef, n, n) + cache = _step_cache(algo, n, optimizer) # Find a suitable initial position - P = P_init - # Iterative procedure + if isnothing(P_init) + P = fill(one(T) / n, n, n) + _flat_gradient!(G, A, B) + else + P = Matrix{T}(P_init) + _gradient!(G, X, A, B, Bt, P, symmetric) + end + # Iterative procedure. The gradient is linear in P, so it can be updated along the + # iterations with a single gradient evaluation (at the search direction Q) per step. + converged = false for _ in 1:max_iter - # Compute the gradient of the quadratic objective at the current point - ∇f = -A * P * B' - A' * P * B - # Compute the search direction Q - Q = graph_matching_search_direction( - algo, ∇f; optimizer, regularization, max_iter_sinkhorn + α, dist_PQ = _frank_wolfe_step!( + algo, + P, + G, + GQ, + X, + A, + B, + Bt, + symmetric, + cache; + regularization, + max_iter_sinkhorn, + tol_sinkhorn, ) - # Compute the step size α - α = graph_matching_step_size(A, B, P, Q) - # Update P - P_new = α * P + (1 - α) * Q - # Check convergence - converged = norm(∇f) * (norm(P - P_new)) < tol - P = P_new + G .= α .* G .+ (1 - α) .* GQ + # Check convergence, knowing that ‖P_new - P‖ = (1 - α) ‖P - Q‖ + converged = (1 - α) * dist_PQ / sqrt(n) < tol converged && break end - # Project onto the set of permutation matrices - P = min_cost_assignment(-P; integer=true, optimizer=optimizer) - # Compute the error - dist = norm(A * P - P * B) + # Project onto the set of permutation matrices: the assignment LP has integral vertices, + # and the integer program is only solved if the optimizer returned a non-vertex solution. + σ = _argmax_permutation(P) + if isnothing(σ) + if isnothing(optimizer) + σ = linear_sum_assignment(P; maximize=true) + else + model = cache isa Model ? cache : _assignment_model(n, optimizer) + σ = _permutation_vector(_max_assignment!(model, P)) + if isnothing(σ) + P_int = min_cost_assignment(-P; integer=true, optimizer=optimizer) + σ = _permutation_vector(P_int) + end + end + end + P = zeros(T, n, n) + for i in 1:n + P[i, σ[i]] = one(T) + end + # Compute the error ‖A P - P B‖ = ‖A - P B Pᵀ‖ + dist = norm(A - B[σ, σ]) return P, dist, converged end diff --git a/src/min_vertex_cover.jl b/src/min_vertex_cover.jl index a92a47e..cf665ac 100644 --- a/src/min_vertex_cover.jl +++ b/src/min_vertex_cover.jl @@ -61,7 +61,7 @@ function min_vertex_cover(g::AbstractGraph; integer::Bool=true, optimizer=HiGHS. @assert termination_status(model) == OPTIMAL cover_values = Vector(value.(model[:cover])) - cover_vertices = findall(==(1), cover_values) + cover_vertices = findall(>(0.5), cover_values) return cover_vertices end diff --git a/test/assignment.jl b/test/assignment.jl index b88b534..acd0e46 100644 --- a/test/assignment.jl +++ b/test/assignment.jl @@ -8,3 +8,15 @@ P = min_cost_assignment(-A; integer=true) @test GraphsOptim.is_doubly_stochastic(min_cost_assignment(rand(10, 10); integer=false)) @test GraphsOptim.is_permutation_matrix(min_cost_assignment(rand(10, 10); integer=true)) + +@testset "Linear sum assignment" begin + for n in (1, 5, 20), + C in (rand(n, n), Float64.(rand(0:2, n, n))), + maximize in (false, true) + + σ = GraphsOptim.linear_sum_assignment(C; maximize) + @test isperm(σ) + P = min_cost_assignment(maximize ? -C : C; integer=true) + @test sum(C[i, σ[i]] for i in 1:n) ≈ sum(C .* P) + end +end diff --git a/test/graph_matching.jl b/test/graph_matching.jl index 85c8f25..cd01b56 100644 --- a/test/graph_matching.jl +++ b/test/graph_matching.jl @@ -1,8 +1,21 @@ using Graphs using GraphsOptim +using HiGHS using LinearAlgebra +using SparseArrays using Test +@testset "Step size" begin + n = 6 + A, B = rand(n, n), rand(n, n) + P = GraphsOptim.flat_doubly_stochastic(n) + Q = Matrix(1.0I, n, n)[[2, 4, 6, 1, 3, 5], :] + f(α) = tr(A * (α * P + (1 - α) * Q) * B' * (α * P + (1 - α) * Q)') + α = GraphsOptim.graph_matching_step_size(A, B, P, Q) + @test 0 <= α <= 1 + @test f(α) >= maximum(f, range(0, 1; length=1001)) - 1e-10 +end + @testset verbose = true "FAQ" begin @test GraphsOptim.is_permutation_matrix( first(graph_matching(FAQ(), rand(5, 5), rand(5, 5))) @@ -38,3 +51,15 @@ end @test converged end end + +@testset "Sparse, dense and JuMP assignment" begin + A, B = sprand(30, 30, 0.2), sprand(30, 30, 0.2) + for algo in (FAQ(), GOAT()) + P_sparse, dist_sparse, _ = graph_matching(algo, A, B) + P_dense, dist_dense, _ = graph_matching(algo, Matrix(A), Matrix(B)) + @test P_sparse == P_dense + @test dist_sparse ≈ dist_dense ≈ norm(A * P_dense - P_dense * B) + end + P_jump, _, _ = graph_matching(FAQ(), A, B; optimizer=HiGHS.Optimizer) + @test GraphsOptim.is_permutation_matrix(P_jump) +end diff --git a/test/min_vertex_cover.jl b/test/min_vertex_cover.jl index 1a5b77e..bae3528 100644 --- a/test/min_vertex_cover.jl +++ b/test/min_vertex_cover.jl @@ -16,3 +16,12 @@ expected_ans_options = [Set([1, 2]), Set([2, 3])] graph = SimpleGraph(adj_matrix) @test Set(min_vertex_cover(graph)) in expected_ans_options + +# Covers must be valid even when the solver returns values slightly below 1 +using Graphs: erdos_renyi, edges, src, dst, nv +for seed in 1:5 + g = erdos_renyi(50, 0.1; seed) + cover = Set(min_vertex_cover(g)) + @test all(src(e) in cover || dst(e) in cover for e in edges(g)) + @test length(cover) == nv(g) - length(maximum_weight_independent_set(g)) +end