Skip to content

Commit

Permalink
update for 0.6 (#64)
Browse files Browse the repository at this point in the history
* update for 0.6

* drop 0.4 support

* modernize and update benchmarks

* fixes

* update benchmark text
  • Loading branch information
KristofferC authored Feb 7, 2017
1 parent 9d09e91 commit 387cfe5
Show file tree
Hide file tree
Showing 13 changed files with 195 additions and 215 deletions.
1 change: 0 additions & 1 deletion .travis.yml
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,6 @@ os:
- linux
- osx
julia:
- 0.4
- 0.5
- nightly
notifications:
Expand Down
143 changes: 73 additions & 70 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,6 @@
[![Build Status](https://travis-ci.org/JuliaStats/Distances.jl.svg?branch=master)](https://travis-ci.org/JuliaStats/Distances.jl)
[![Coverage Status](https://coveralls.io/repos/JuliaStats/Distances.jl/badge.svg?branch=master&service=github)](https://coveralls.io/github/JuliaStats/Distances.jl?branch=master)

[![Distances](http://pkg.julialang.org/badges/Distances_0.4.svg)](http://pkg.julialang.org/?pkg=Distances&ver=0.4)
[![Distances](http://pkg.julialang.org/badges/Distances_0.5.svg)](http://pkg.julialang.org/?pkg=Distances)

A Julia package for evaluating distances(metrics) between vectors.
Expand Down Expand Up @@ -125,32 +124,32 @@ This type system has practical significance. For example, when computing pairwis

Each distance corresponds to a distance type. The type name and the corresponding mathematical definitions of the distances are listed in the following table.

| type name | convenient syntax | math definition |
| -------------------- | ------------------------ | --------------------|
| Euclidean | euclidean(x, y) | sqrt(sum((x - y) .^ 2)) |
| SqEuclidean | sqeuclidean(x, y) | sum((x - y).^2) |
| Cityblock | cityblock(x, y) | sum(abs(x - y)) |
| Chebyshev | chebyshev(x, y) | max(abs(x - y)) |
| Minkowski | minkowski(x, y, p) | sum(abs(x - y).^p) ^ (1/p) |
| Hamming | hamming(x, y) | sum(x .!= y) |
| Rogers-Tanimoto | rogerstanimoto(x, y) | 2(sum(x&!y) + sum(!x&y)) / (2(sum(x&!y) + sum(!x&y)) + sum(x&y) + sum(!x&!y)) |
| Jaccard | jaccard(x, y) | 1 - sum(min(x, y)) / sum(max(x, y)) |
| CosineDist | cosine_dist(x, y) | 1 - dot(x, y) / (norm(x) * norm(y)) |
| CorrDist | corr_dist(x, y) | cosine_dist(x - mean(x), y - mean(y)) |
| ChiSqDist | chisq_dist(x, y) | sum((x - y).^2 / (x + y)) |
| KLDivergence | kl_divergence(x, y) | sum(p .* log(p ./ q)) |
| RenyiDivergence | renyi_divergence(x, y, k)| log(sum( x .* (x ./ y) .^ (k - 1))) / (k - 1) |
| JSDivergence | js_divergence(x, y) | KL(x, m) / 2 + KL(y, m) / 2 with m = (x + y) / 2 |
| SpanNormDist | spannorm_dist(x, y) | max(x - y) - min(x - y ) |
| BhattacharyyaDist | bhattacharyya(x, y) | -log(sum(sqrt(x .* y) / sqrt(sum(x) * sum(y))) |
| HellingerDist | hellinger(x, y) | sqrt(1 - sum(sqrt(x .* y) / sqrt(sum(x) * sum(y)))) |
| Mahalanobis | mahalanobis(x, y, Q) | sqrt((x - y)' * Q * (x - y)) |
| SqMahalanobis | sqmahalanobis(x, y, Q) | (x - y)' * Q * (x - y) |
| WeightedEuclidean | weuclidean(x, y, w) | sqrt(sum((x - y).^2 .* w)) |
| WeightedSqEuclidean | wsqeuclidean(x, y, w) | sum((x - y).^2 .* w) |
| WeightedCityblock | wcityblock(x, y, w) | sum(abs(x - y) .* w) |
| WeightedMinkowski | wminkowski(x, y, w, p) | sum(abs(x - y).^p .* w) ^ (1/p) |
| WeightedHamming | whamming(x, y, w) | sum((x .!= y) .* w) |
| type name | convenient syntax | math definition |
| -------------------- | -------------------------- | --------------------|
| Euclidean | `euclidean(x, y)` | `sqrt(sum((x - y) .^ 2))` |
| SqEuclidean | `sqeuclidean(x, y)` | `sum((x - y).^2)` |
| Cityblock | `cityblock(x, y)` | `sum(abs(x - y))` |
| Chebyshev | `chebyshev(x, y)` | `max(abs(x - y))` |
| Minkowski | `minkowski(x, y, p)` | `sum(abs(x - y).^p) ^ (1/p)` |
| Hamming | `hamming(x, y)` | `sum(x .!= y)` |
| Rogers-Tanimoto | `rogerstanimoto(x, y)` | `2(sum(x&!y) + sum(!x&y)) / (2(sum(x&!y) + sum(!x&y)) + sum(x&y) + sum(!x&!y))` |
| Jaccard | `jaccard(x, y)` | `1 - sum(min(x, y)) / sum(max(x, y))` |
| CosineDist | `cosine_dist(x, y)` | `1 - dot(x, y) / (norm(x) * norm(y))` |
| CorrDist | `corr_dist(x, y)` | `cosine_dist(x - mean(x), y - mean(y))` |
| ChiSqDist | `chisq_dist(x, y)` | `sum((x - y).^2 / (x + y))` |
| KLDivergence | `kl_divergence(x, y)` | `sum(p .* log(p ./ q))` |
| RenyiDivergence | `renyi_divergence(x, y, k)`| `log(sum( x .* (x ./ y) .^ (k - 1))) / (k - 1)` |
| JSDivergence | `js_divergence(x, y)` | `KL(x, m) / 2 + KL(y, m) / 2 with m = (x + y) / 2` |
| SpanNormDist | `spannorm_dist(x, y)` | `max(x - y) - min(x - y )` |
| BhattacharyyaDist | `bhattacharyya(x, y)` | `-log(sum(sqrt(x .* y) / sqrt(sum(x) * sum(y)))` |
| HellingerDist | `hellinger(x, y) ` | `sqrt(1 - sum(sqrt(x .* y) / sqrt(sum(x) * sum(y))))` |
| Mahalanobis | `mahalanobis(x, y, Q)` | `sqrt((x - y)' * Q * (x - y))` |
| SqMahalanobis | `sqmahalanobis(x, y, Q)` | ` (x - y)' * Q * (x - y)` |
| WeightedEuclidean | `weuclidean(x, y, w)` | `sqrt(sum((x - y).^2 .* w))` |
| WeightedSqEuclidean | `wsqeuclidean(x, y, w)` | `sum((x - y).^2 .* w)` |
| WeightedCityblock | `wcityblock(x, y, w)` | `sum(abs(x - y) .* w)` |
| WeightedMinkowski | `wminkowski(x, y, w, p)` | `sum(abs(x - y).^p .* w) ^ (1/p)` |
| WeightedHamming | `whamming(x, y, w)` | `sum((x .!= y) .* w)` |

**Note:** The formulas above are using *Julia*'s functions. These formulas are mainly for conveying the math concepts in a concise way. The actual implementation may use a faster way.

Expand Down Expand Up @@ -186,58 +185,62 @@ julia> pairwise(Euclidean(1e-12), x, x)

The implementation has been carefully optimized based on benchmarks. The Julia scripts ``test/bench_colwise.jl`` and ``test/bench_pairwise.jl`` run the benchmarks on a variety of distances, respectively under column-wise and pairwise settings.

Here are the benchmarks that I obtained on Mac OS X 10.8 with Intel Core i7 2.6 GHz.
Here are benchmarks obtained on Linux with Intel Core i7-4770K 3.5 GHz.

#### Column-wise benchmark

The table below compares the performance (measured in terms of average elapsed time of each iteration) of a straightforward loop implementation and an optimized implementation provided in *Distances.jl*. The task in each iteration is to compute a specific distance between corresponding columns in two ``200-by-10000`` matrices.

| distance | loop | colwise | gain |
|------------ | --------| ------------| -----------|
| SqEuclidean | 0.046962 | 0.002782 | 16.8782 |
| Euclidean | 0.046667 | 0.0029 | 16.0937 |
| Cityblock | 0.046619 | 0.0031 | 15.039 |
| Chebyshev | 0.053578 | 0.010856 | 4.9356 |
| Minkowski | 0.061804 | 0.02357 | 2.6221 |
| Hamming | 0.044047 | 0.00219 | 20.1131 |
| CosineDist | 0.04496 | 0.002855 | 15.7457 |
| CorrDist | 0.080828 | 0.029708 | 2.7207 |
| ChiSqDist | 0.051009 | 0.008088 | 6.307 |
| KLDivergence | 0.079598 | 0.035353 | 2.2515 |
| JSDivergence | 0.545789 | 0.493362 | 1.1063 |
| WeightedSqEuclidean | 0.046182 | 0.003219 | 14.3477 |
| WeightedEuclidean | 0.046831 | 0.004122 | 11.3603 |
| WeightedCityblock | 0.046457 | 0.003636 | 12.7781 |
| WeightedMinkowski | 0.062532 | 0.020486 | 3.0524 |
| WeightedHamming | 0.046217 | 0.002269 | 20.3667 |
| SqMahalanobis | 0.150364 | 0.042335 | 3.5518 |
| Mahalanobis | 0.159638 | 0.041071 | 3.8869 |

We can see that using ``colwise`` instead of a simple loop yields considerable gain (2x - 9x), especially when the internal computation of each distance is simple. Nonetheless, when the computaton of a single distance is heavy enough (e.g. *Minkowski* and *JSDivergence*), the gain is not as significant.
| distance | loop | colwise | gain |
|----------- | -------| ----------| -------|
| SqEuclidean | 0.012308s | 0.003860s | 3.1884 |
| Euclidean | 0.012484s | 0.003995s | 3.1246 |
| Cityblock | 0.012463s | 0.003927s | 3.1735 |
| Chebyshev | 0.014897s | 0.005898s | 2.5258 |
| Minkowski | 0.028154s | 0.017812s | 1.5806 |
| Hamming | 0.012200s | 0.003896s | 3.1317 |
| CosineDist | 0.013816s | 0.004670s | 2.9583 |
| CorrDist | 0.023349s | 0.016626s | 1.4044 |
| ChiSqDist | 0.015375s | 0.004788s | 3.2109 |
| KLDivergence | 0.044360s | 0.036123s | 1.2280 |
| JSDivergence | 0.098587s | 0.085595s | 1.1518 |
| BhattacharyyaDist | 0.023103s | 0.013002s | 1.7769 |
| HellingerDist | 0.023329s | 0.012555s | 1.8581 |
| WeightedSqEuclidean | 0.012136s | 0.003758s | 3.2296 |
| WeightedEuclidean | 0.012307s | 0.003789s | 3.2482 |
| WeightedCityblock | 0.012287s | 0.003923s | 3.1321 |
| WeightedMinkowski | 0.029895s | 0.018471s | 1.6185 |
| WeightedHamming | 0.013427s | 0.004082s | 3.2896 |
| SqMahalanobis | 0.121636s | 0.019370s | 6.2796 |
| Mahalanobis | 0.117871s | 0.019939s | 5.9117 |

We can see that using ``colwise`` instead of a simple loop yields considerable gain (2x - 6x), especially when the internal computation of each distance is simple. Nonetheless, when the computaton of a single distance is heavy enough (e.g. *Minkowski* and *JSDivergence*), the gain is not as significant.

#### Pairwise benchmark

The table below compares the performance (measured in terms of average elapsed time of each iteration) of a straightforward loop implementation and an optimized implementation provided in *Distances.jl*. The task in each iteration is to compute a specific distance in a pairwise manner between columns in a ``100-by-200`` and ``100-by-250`` matrices, which will result in a ``200-by-250`` distance matrix.

| distance | loop | pairwise | gain |
|------------ | --------| ------------| -----------|
| SqEuclidean | 0.119961 | 0.00037 | **324.6457** |
| Euclidean | 0.122645 | 0.000678 | **180.9180** |
| Cityblock | 0.116956 | 0.007997 | 14.6251 |
| Chebyshev | 0.137985 | 0.028489 | 4.8434 |
| Minkowski | 0.170101 | 0.059991 | 2.8354 |
| Hamming | 0.110742 | 0.004781 | 23.1627 |
| CosineDist | 0.110913 | 0.000514 | **215.8028** |
| CorrDist | 0.1992 | 0.000808 | 246.4574 |
| ChiSqDist | 0.124782 | 0.020781 | 6.0046 |
| KLDivergence | 0.1994 | 0.088366 | 2.2565 |
| JSDivergence | 1.35502 | 1.215785 | 1.1145 |
| WeightedSqEuclidean | 0.119797 | 0.000444 | **269.531** |
| WeightedEuclidean | 0.126304 | 0.000712 | **177.5122** |
| WeightedCityblock | 0.117185 | 0.011475 | 10.2122 |
| WeightedMinkowski | 0.172614 | 0.061693 | 2.7979 |
| WeightedHamming | 0.112525 | 0.005072 | 22.1871 |
| SqMahalanobis | 0.377342 | 0.000577 | **653.9759** |
| Mahalanobis | 0.373796 | 0.002359 | **158.4337** |
| distance | loop | pairwise | gain |
|----------- | -------| ----------| -------|
| SqEuclidean | 0.032179s | 0.000170s | **189.7468** |
| Euclidean | 0.031646s | 0.000326s | **97.1773** |
| Cityblock | 0.031594s | 0.002771s | 11.4032 |
| Chebyshev | 0.036732s | 0.011575s | 3.1735 |
| Minkowski | 0.073685s | 0.047725s | 1.5440 |
| Hamming | 0.030016s | 0.002539s | 11.8236 |
| CosineDist | 0.035426s | 0.000235s | **150.8504** |
| CorrDist | 0.061430s | 0.000341s | **180.1693** |
| ChiSqDist | 0.037702s | 0.011709s | 3.2199 |
| KLDivergence | 0.119043s | 0.086861s | 1.3705 |
| JSDivergence | 0.255449s | 0.227079s | 1.1249 |
| BhattacharyyaDist | 0.059165s | 0.033330s | 1.7751 |
| HellingerDist | 0.056953s | 0.031163s | 1.8276 |
| WeightedSqEuclidean | 0.031781s | 0.000218s | **145.9820** |
| WeightedEuclidean | 0.031365s | 0.000410s | **76.4517** |
| WeightedCityblock | 0.031239s | 0.003242s | 9.6360 |
| WeightedMinkowski | 0.077039s | 0.049319s | 1.5621 |
| WeightedHamming | 0.032584s | 0.005673s | 5.7442 |
| SqMahalanobis | 0.280485s | 0.000297s | **943.6018** |
| Mahalanobis | 0.295715s | 0.000498s | **593.6096** |

For distances of which a major part of the computation is a quadratic form (e.g. *Euclidean*, *CosineDist*, *Mahalanobis*), the performance can be drastically improved by restructuring the computation and delegating the core part to ``GEMM`` in *BLAS*. The use of this strategy can easily lead to 100x performance gain over simple loops (see the highlighted part of the table above).
3 changes: 1 addition & 2 deletions REQUIRE
Original file line number Diff line number Diff line change
@@ -1,2 +1 @@
julia 0.4
Compat 0.8.4
julia 0.5
2 changes: 0 additions & 2 deletions src/Distances.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,6 @@ __precompile__()

module Distances

import Compat.view

export
# generic types/functions
PreMetric,
Expand Down
6 changes: 3 additions & 3 deletions src/common.jl
Original file line number Diff line number Diff line change
Expand Up @@ -101,7 +101,7 @@ end
function sumsq_percol{T}(a::AbstractMatrix{T})
m = size(a, 1)
n = size(a, 2)
r = Array(T, n)
r = Vector{T}(n)
for j = 1:n
aj = view(a, :, j)
r[j] = dot(aj, aj)
Expand All @@ -113,7 +113,7 @@ function wsumsq_percol{T1, T2}(w::AbstractArray{T1}, a::AbstractMatrix{T2})
m = size(a, 1)
n = size(a, 2)
T = typeof(one(T1)*one(T2))
r = Array(T, n)
r = Vector{T}(n)
for j = 1:n
aj = view(a, :, j)
s = zero(T)
Expand All @@ -138,4 +138,4 @@ function dot_percol!(r::AbstractArray, a::AbstractMatrix, b::AbstractMatrix)
return r
end

dot_percol(a::AbstractMatrix, b::AbstractMatrix) = dot_percol!(Array(Float64, size(a,2)), a, b)
dot_percol(a::AbstractMatrix, b::AbstractMatrix) = dot_percol!(Vector{Float64}(size(a,2)), a, b)
16 changes: 5 additions & 11 deletions src/generic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -62,19 +62,19 @@ end

function colwise(metric::PreMetric, a::AbstractMatrix, b::AbstractMatrix)
n = get_common_ncols(a, b)
r = Array(result_type(metric, a, b), n)
r = Vector{result_type(metric, a, b)}(n)
colwise!(r, metric, a, b)
end

function colwise(metric::PreMetric, a::AbstractVector, b::AbstractMatrix)
n = size(b, 2)
r = Array(result_type(metric, a, b), n)
r = Vector{result_type(metric, a, b)}(n)
colwise!(r, metric, a, b)
end

function colwise(metric::PreMetric, a::AbstractMatrix, b::AbstractVector)
n = size(a, 2)
r = Array(result_type(metric, a, b), n)
r = Vector{result_type(metric, a, b)}(n)
colwise!(r, metric, a, b)
end

Expand Down Expand Up @@ -117,18 +117,12 @@ end
function pairwise(metric::PreMetric, a::AbstractMatrix, b::AbstractMatrix)
m = size(a, 2)
n = size(b, 2)
r = Array(result_type(metric, a, b), (m, n))
r = Matrix{result_type(metric, a, b)}(m, n)
pairwise!(r, metric, a, b)
end

function pairwise(metric::PreMetric, a::AbstractMatrix)
n = size(a, 2)
r = Array(result_type(metric, a, a), (n, n))
pairwise!(r, metric, a)
end

function pairwise(metric::SemiMetric, a::AbstractMatrix)
n = size(a, 2)
r = Array(result_type(metric, a, a), (n, n))
r = Matrix{result_type(metric, a, a)}(n, n)
pairwise!(r, metric, a)
end
4 changes: 2 additions & 2 deletions src/metrics.jl
Original file line number Diff line number Diff line change
Expand Up @@ -341,8 +341,8 @@ rogerstanimoto{T <: Bool}(a::AbstractArray{T}, b::AbstractArray{T}) = evaluate(R
# SqEuclidean
function pairwise!(r::AbstractMatrix, dist::SqEuclidean, a::AbstractMatrix, b::AbstractMatrix)
At_mul_B!(r, a, b)
sa2 = sumabs2(a, 1)
sb2 = sumabs2(b, 1)
sa2 = sum(abs2, a, 1)
sb2 = sum(abs2, b, 1)
threshT = convert(eltype(r), dist.thresh)
if threshT <= 0
# If there's no chance of triggering the threshold, we can use @simd
Expand Down
1 change: 0 additions & 1 deletion test/REQUIRE

This file was deleted.

87 changes: 41 additions & 46 deletions test/bench_colwise.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,34 +2,26 @@
# Benchmark on column-wise distance evaluation

using Distances
using BenchmarkTools

macro bench_colwise_dist(repeat, dist, x, y)
quote
println("bench ", typeof($dist))
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 1.0

# warming up
r1 = evaluate($dist, ($x)[:,1], ($y)[:,1])
colwise($dist, $x, $y)
function bench_colwise_distance(dist, x, y)
r1 = evaluate(dist, x[:,1], y[:,1])

# timing

t0 = @elapsed for k = 1 : $repeat
n = size($x, 2)
r = Array(typeof(r1), n)
for j = 1 : n
r[j] = evaluate($dist, ($x)[:,j], ($y)[:,j])
end
end
@printf " loop: t = %9.6fs\n" (t0 / $repeat)

t1 = @elapsed for k = 1 : $repeat
r = colwise($dist, $x, $y)
# timing
t0 = @belapsed begin
n = size(x, 2)
r = Vector{typeof($r1)}(n)
for j = 1:n
r[j] = evaluate($dist, $(x)[:, j], $(y)[:, j])
end
@printf " colwise: t = %9.6fs | gain = %7.4fx\n" (t1 / $repeat) (t0 / t1)
println()
end
end

t1 = @belapsed colwise($dist, $x, $y)
print("| ", typeof(dist).name.name, " |")
@printf("%9.6fs | %9.6fs | %7.4f |\n", t0, t1, (t0 / t1))
end

m = 200
n = 10000
Expand All @@ -41,27 +33,30 @@ w = rand(m)
Q = rand(m, m)
Q = Q' * Q

@bench_colwise_dist 20 SqEuclidean() x y
@bench_colwise_dist 20 Euclidean() x y
@bench_colwise_dist 20 Cityblock() x y
@bench_colwise_dist 20 Chebyshev() x y
@bench_colwise_dist 5 Minkowski(3.0) x y
@bench_colwise_dist 20 Hamming() x y

@bench_colwise_dist 20 CosineDist() x y
@bench_colwise_dist 10 CorrDist() x y
@bench_colwise_dist 20 ChiSqDist() x y
@bench_colwise_dist 10 KLDivergence() x y
@bench_colwise_dist 5 JSDivergence() x y

@bench_colwise_dist 10 BhattacharyyaDist() x y
@bench_colwise_dist 10 HellingerDist() x y

@bench_colwise_dist 20 WeightedSqEuclidean(w) x y
@bench_colwise_dist 20 WeightedEuclidean(w) x y
@bench_colwise_dist 20 WeightedCityblock(w) x y
@bench_colwise_dist 5 WeightedMinkowski(w, 3.0) x y
@bench_colwise_dist 20 WeightedHamming(w) x y

@bench_colwise_dist 10 SqMahalanobis(Q) x y
@bench_colwise_dist 10 Mahalanobis(Q) x y
println("| distance | loop | colwise | gain |")
println("|----------- | -------| ----------| -------|")

bench_colwise_distance(SqEuclidean(), x, y)
bench_colwise_distance(Euclidean(), x, y)
bench_colwise_distance(Cityblock(), x, y)
bench_colwise_distance(Chebyshev(), x, y)
bench_colwise_distance(Minkowski(3.0), x, y)
bench_colwise_distance(Hamming(), x, y)

bench_colwise_distance(CosineDist(), x, y)
bench_colwise_distance(CorrDist(), x, y)
bench_colwise_distance(ChiSqDist(), x, y)
bench_colwise_distance(KLDivergence(), x, y)
bench_colwise_distance(JSDivergence(), x, y)

bench_colwise_distance(BhattacharyyaDist(), x, y)
bench_colwise_distance(HellingerDist(), x, y)

bench_colwise_distance(WeightedSqEuclidean(w), x, y)
bench_colwise_distance(WeightedEuclidean(w), x, y)
bench_colwise_distance(WeightedCityblock(w), x, y)
bench_colwise_distance(WeightedMinkowski(w, 3.0), x, y)
bench_colwise_distance(WeightedHamming(w), x, y)

bench_colwise_distance(SqMahalanobis(Q), x, y)
bench_colwise_distance(Mahalanobis(Q), x, y)
Loading

0 comments on commit 387cfe5

Please sign in to comment.