# Fast hamming distance function in Julia that returns a distance matrix

I would like to use Julia to compute the hamming distance on a very large dataset. I need to get back a distance matrix between rows in order to run further analysis on this matrix.

For my purposes, it is useful that the data are stored in a DataFrame type.

using DataFrames


The data looks something like this

a = [1 0 1 0 ; 1 1 1 1; 0 0 0 0; 0 0 0 0 ; 0 0 0 1]

df = convert(DataFrame, a);

nrows = size(df, 1)
ncols = size(df, 2)


I made a function in Julia

function hamjulia(df)

nrows = size(df, 1)
ncols = size(df, 2)
m, n = nrows, nrows
A = fill(0, (m, n))

for i in 1:nrows
for k in 1:nrows

v = 0
for j in 1:ncols
if df[i,j] != df[k,j]
v += 1
end
end
A[i,k] = v
end
end
return A
end

p = hamjulia(df)
p


My issue is that this code is super slow compared to some R packages. For instance, when I compared this function to the rdist package, rdist(df, metric = 'hamming'), R is faster.

How could I make this code really efficient? Especially that I would need to run it on a very large dataframe.

Thanks.

I have three areas of suggestions below. Here's a modified version of your code that implements these suggestions:

# For memory safety, restrict type to Matrix due to the use of @inbounds below.
function ham(df::Matrix)
nrows, ncols = size(df)
A = fill(0, ncols, ncols)
for j in 1:ncols, i in j+1:ncols
v = 0
for k in 1:nrows
v += @inbounds (df[k, i] != df[k, j])
end
A[i, j] = A[j, i] = v
end
return A
end

ham(df::DataFrame) = ham(convert(Matrix, df))


1. Write code for column-major arrays.

Julia is column-major (and packages like DataFrames.jl conform with this by storing data in columns), so when performance is important, you should iterate through the rows of a column. For example, when iterating through an array A[x, y, z], the leftmost index x should be changing the most rapidly. Looking at the innermost loop in your code,

for j in 1:ncols
if df[i,j] != df[k,j]
v += 1
end
end


notice that it is iterating through the columns of a row instead (i.e., j is changing most rapidly in df[i,j]). To get a free performance boost, rewrite the code to exchange the roles of columns and rows. The data will also need to be transposed. (Note that R is also column-major.)

Julia documentation: memory order.

2. In performance-critical loops, annotate @inbounds where possible.

In the modified code I gave, the array access in the innermost loop is annotated with @inbounds to disable bounds checking. This gives a speedup in large part because it allows compiler to automatically vectorize the loop.

Only use @inbounds if you are certain that out-of-bounds access is impossible; in our case, we cannot guarantee this for custom types, so we restrict the input type to a Matrix (an alias for Array{T,2} where T) to ensure that df[k, i] and df[k, j] are always in bounds. Like so:

function ham(df::Matrix)
[...]
for k in 1:nrows
v += @inbounds (df[k, i] != df[k, j])
end
[...]
end


If @inbounds is not used, I would not restrict the input type and would instead leave it generic.

Julia documentation: performance annotations.

3. Beware of type instability.

For performant code, functions should be type stable, meaning that given an input type, the output type can be inferred during compilation. Similarly, variables used within a function should be type stable.

As noted in its documentation, DataFrame columns are not type stable! That is, the type of a DataFrame does not contain information about the types of its columns; in fact, the DataFrame type is not parametric, so the type contains no information other than that it is a DataFrame. This can drastically degrade performance.

The @code_warntype macro is very helpful for diagnosing this.

In our case, assuming that all columns are of the same type, a simple workaround is to copy the columns into a Matrix before calling the function. (The time this takes is negligible compared to the rest of the function.) To do so, we can specialize ham for the DataFrame type:

ham(df::DataFrame) = ham(convert(Matrix, df))


Note that the compiler specializes code for types at function boundaries, so it is important that this conversion takes place outside the original ham function. In particular, don't put the conversion into the original function like so:

# Don't do this!
function ham(x)
df = convert(Matrix, x)
[rest of ham...]
end


because the type of x here might not contain sufficient information to infer the type of df. Such is the case for DataFrames.

Julia documentation: function barriers.

Simple benchmark

After implementing these suggestions, the speedup on calculating the distance matrix for 128 vectors of 256 integers each is by about 1000x: (note that the times are in microseconds vs milliseconds)

using BenchmarkTools
using DataFrames
M = rand(0:1, 256, 128);
df = DataFrame(M);
tdf = DataFrame(permutedims(M)); # transpose rows and columns

julia> @benchmark ham($df) BenchmarkTools.Trial: memory estimate: 387.33 KiB allocs estimate: 133 -------------- minimum time: 493.841 μs (0.00% GC) median time: 537.689 μs (0.00% GC) mean time: 610.732 μs (2.74% GC) maximum time: 6.215 ms (81.69% GC) -------------- samples: 8176 evals/sample: 1  julia> @benchmark hamjulia($tdf)
BenchmarkTools.Trial:
memory estimate:  128.08 KiB
allocs estimate:  2
--------------
minimum time:     538.585 ms (0.00% GC)
median time:      546.390 ms (0.00% GC)
mean time:        551.220 ms (0.00% GC)
maximum time:     567.238 ms (0.00% GC)
--------------
samples:          10
evals/sample:     1


The additional memory allocations in ham are due to the conversion of the DataFrame to a Matrix.

• Thanks a lot for this! I will run some tests with it. I noticed that I don't need the DataFrame so much. What about only filling half of the distance matrix? Do you know how stop when we reached the diagonal? (the distance matrix is filled twice here) – giac Mar 17 at 9:48
• No problem! The sample code already stops at the diagonal: notice that the outer loop goes for j in 1:ncols, i in j+1:ncols, where i depends on j and so excludes the diagonal and one half of the matrix. Hence the A[i, j] = A[j, i] = v sets each non-diagonal element exactly once. – Vincent Yu Mar 17 at 10:04
• Oh sorry I missed that. But your function still returns a full matrix right? Wouldn't save memory to only return the half matrix? – giac Mar 17 at 10:33
• Yes, it returns a full matrix. – Vincent Yu Mar 17 at 10:34
• Sorry to bother you with this, but it return the full matrix but only compute the half matrix. Is that correct? Thanks – giac Mar 17 at 10:37

As a cheap first try I suggest Distances.jl with (I think)

pairwise(Hamming(), df, dims=1)