How to parallelize computation of pairwise distance matrix?

Question

My problem is roughly as follows. Given a numerical matrix X, where each row is an item. I want to find each row's nearest neighbor in terms of L2 distance in all rows except itself. I tried reading the official documentation but was still a little confused about how to achieve this. Could someone give me some hint?

My code is as follows

function l2_dist(v1, v2)
    return sqrt(sum((v1 - v2) .^ 2))
end

function main(Mat, dist_fun)
    n = size(Mat, 1)

    Dist = SharedArray{Float64}(n) #[Inf for i in 1:n]
    Id = SharedArray{Int64}(n) #[-1 for i in 1:n]
    @parallel for i = 1:n
        Dist[i] = Inf
        Id[i] = 0
    end
    Threads.@threads for i in 1:n
        for j in 1:n
            if i != j
                println(i, j)
                dist_temp = dist_fun(Mat[i, :], Mat[j, :])
                if dist_temp < Dist[i]
                    println("Dist updated!")
                    Dist[i] = dist_temp
                    Id[i] = j
                end
            end
        end
    end
    return Dict("Dist" => Dist, "Id" => Id)
end

n = 4000
p = 30

X = [rand() for i in 1:n, j in 1:p];


main(X[1:30, :], l2_dist)
@time N = main(X, l2_dist)

I'm trying to distributed all the i's (i.e. calculating each row minimum) over different cores. But the version above apparently isn't working correctly. It is even slower than the sequential version. Can someone point me to the right direction? Thanks.

Answer 1

Maybe you're doing something in addition to what you have written down, but, at this point from what I can see, you aren't actually doing any computations in parallel. Julia requires you to tell it how many processors (or threads) you would like it to have access to. You can do this through either

• Starting Julia with multiple processors julia -p # (where # is the number of processors you want Julia to have access to)
• Once you have started a Julia "session" you can call the addprocs function to add additional processors.
• To have more than 1 thread, you need to run command export JULIA_NUM_THREADS = #. I don't know very much about threading, so I will be sticking with the @parallel macro. I suggest reading documentation for more details on threading -- Maybe @Chris Rackauckas could expand a little more on the difference.

A few comments below about my code and on your code:

• I'm on version 0.6.1-pre.0. I don't think I'm doing anything 0.6 specific, but this is a heads up just in case.
• I'm going to use the Distances.jl package when computing the distances between vectors. I think it is a good habit to farm out as many of my computations to well-written and well-maintained packages as possible.
• Rather than compute the distance between rows, I'm going to compute the distance between columns. This is because Julia is a column-major language, so this will increase the number of cache hits and give a little extra speed. You can obviously get the row-wise results you want by just transposing the input.
• Unless you expect to have that many memory allocations then that many allocations are a sign that something in your code is inefficient. It is often a type stability problem. I don't know if that was the case in your code before, but that doesn't seem to be an issue in the current version (it wasn't immediately clear to me why you were having so many allocations).

Code is below

# Make sure all processors have access to Distances package
@everywhere using Distances

# Create a random matrix
nrow = 30
ncol = 4000

# Seed creation of random matrix so it is always same matrix
srand(42)
X = rand(nrow, ncol)

function main(X::AbstractMatrix{Float64}, M::Distances.Metric)
    # Get size of the matrix
    nrow, ncol = size(X)

    # Create `SharedArray` to store output
    ind_vec = SharedArray{Int}(ncol)
    dist_vec = SharedArray{Float64}(ncol)

    # Compute the distance between columns
    @sync @parallel for i in 1:ncol
        # Initialize various temporary variables
        min_dist_i = Inf
        min_ind_i = -1
        X_i = view(X, :, i)

        # Check distance against all other columns
        for j in 1:ncol
            # Skip comparison with itself
            if i==j
                continue
            end

            # Tell us who is doing the work 
            # (can uncomment if you want to verify stuff)
            # println("Column $i compared with Column $j by worker $(myid())")

            # Evaluate the new distance... 
            # If it is less then replace it, otherwise proceed
            dist_temp = evaluate(M, X_i, view(X, :, j))
            if dist_temp < min_dist_i
                min_dist_i = dist_temp
                min_ind_i = j
            end
        end

        # Which column is minimum distance from column i
        dist_vec[i] = min_dist_i
        ind_vec[i] = min_ind_i
    end

    return dist_vec, ind_vec
end

# Using Euclidean metric
metric = Euclidean()
inds, dist = main(X, metric)

@time main(X, metric);
@show dist[[1, 5, 25]], inds[[1, 5, 25]]

You can run the code with

• 1 processor julia testfile.jl

  % julia testfile.jl
    0.640365 seconds (16.00 M allocations: 732.495 MiB, 3.70% gc time)
    (dist[[1, 5, 25]], inds[[1, 5, 25]]) = ([2541, 2459, 1602], [1.40892, 1.38206, 1.32184])
  

• n processors (in this case 4) julia -p n testfile.jl

  % julia -p 4 testfile.jl
    0.201523 seconds (2.10 k allocations: 99.107 KiB)
    (dist[[1, 5, 25]], inds[[1, 5, 25]]) = ([2541, 2459, 1602], [1.40892, 1.38206, 1.32184])