Parallelizing step3

[matbesancon]

We're using Hungarian.jl in the https://github.com/ZIB-IOL/FrankWolfe.jl package to solve LPs over the Birkhoff polytope.

It scales well in terms of memory footprint but badly in terms of time, mostly because of step3!, where most of the time seems to be spent in this block:

    @inbounds for j in colUncoveredIdx, i in rowUncoveredIdx
        cost = costMat[i, j] + Δcol[j] + Δrow[i]
        if cost <= h
            if cost != h
                h = cost
                empty!(minLocations)
            end
            push!(minLocations, (i, j))
        end
    end

this would be non-trivial to parallelize because of the shared bound h and array minLocations but could be worth it.

[Gnimuc]

It scales well in terms of memory footprint but badly in terms of time, mostly because of step3!, where most of the time seems to be spent in this block:

I think the hot-pot depends on the problem size. Yes, for small matrices, step3 is the bottleneck, but I'm not sure how much performance gain we can benefit from multi-threading.

[Gnimuc]

Looks like I tried something similar years ago:

Hungarian.jl/src/Munkres.jl

Line 343 in a60d303

@threads for j in colUncoveredIdx

[matbesancon]

Do you get other crititcal points when matrix sizes grow larger? This was for a 1000x1000 dense matrix

[Gnimuc]

Hungarian.jl/src/Munkres.jl

Lines 91 to 108 in fe5c402

	for ii in CartesianIndices(size(costMat))
	# "consider a zero Z of the matrix;"
	r, c = ii.I
	cost = costMat[r, c] + Δrow[r] + Δcol[c]
	if cost == 0
	Zs[r, c] = Z
	# "if there is no starred zero in its row and none in its column, star Z.
	# repeat, considering each zero in the matrix in turn;"
	if !colSTAR[c] && !rowSTAR[r]
	Zs[r, c] = STAR
	rowSTAR[r] = true
	colSTAR[c] = true
	row2colSTAR[r] = c
	# "then cover every column containing a starred zero."
	colCovered[c] = true
	end
	end
	end

This marking step is critical for large matrices like 4000x4000 or 8000x8000. You can find some benchmark reports in #15.

[Gnimuc]

What's the problem size of your use case?

[matbesancon]

the problem go up to npixels * npixels for image problems so as big as you want, 150528 * 150528 in some cases but we couldn't run Hungarian.jl for these sizes

[Gnimuc]

That's really large. Hungarian is an O(n^3) algorithm and may not suite for this kinda problem.

Would you like to make a PR with above-mentioned multi-threading support? I'm happy to merge if the benchmark reports are green.

BTW, this is released under public domain, so you could also make the code tailored to your use case and embed it in your package.

[matbesancon]

That's really large. Hungarian is an O(n^3) algorithm and may not suite for this kinda problem.

Yes this is the kind of extreme limit we'd like to solve, it's okay to have the method work for smaller sizes of course.

I'll probably try to reimplement the multithreaded version and compare the two + a LP solver