Performance Bottle Necks

[JonasIsensee]

Most of the time both during searches and graph creation
is spent within distance computation and insert! & pop_nearest!
of the NeighborSets.

To reach the speed of the C++ implementation ( https://github.com/nmslib/hnsw )
both parts need to be sped up.
AFAICT the algorithm does the same as the C++ implementation so
my version shouldn't be doing more unnecessary iterations which
would have been the simplest explanation.

One way to speed up the process should be to use PriorityQueues
instead of my NeighborSets as the C++ uses these. All attempts with PriorityQueue have so far not been faster though. ( I used the PriorityQueue from DataStructures.jl )

The second point is the distance computation. I have no idea if it is possible or desirable to optimize this further. ATM we use the Distances.jl library.
The relevant code in C++ is here: https://github.com/nmslib/hnsw/blob/master/hnswlib/space_l2.h

EDIT: I investigated the @simd optimizations:
in hnsw the runtime is >2× longer for dim % 16 != 0 && dim % 4 != 0 as they implement simd
optimization for those two cases only.
In Julia I do not see this behaviour at all. Is the @simd macro so powerful that it that it works with other dimensions as well? (by adding 0s at the end?)

EDIT: PriorityQueue from DataStructures.jl is very old code in Julian terms (according to the git blame). Maybe it could be optimized now?

[JonasIsensee]

@zgornel
You mentioned that you might be able to help with some of the optimization.
I'm somewhat stuck.

If you have the time, please have a look at my work and if you see anything
don't hesitate to point it out.
If you need additional info on the algorithm, please say so.

[zgornel]

@JonasIsensee I will (been tracking silently the progress). Will put some effort next week (still reading the paper). Anyway, nice work! By the way, the performance issues you mention are for the training phase ?

[JonasIsensee]

Great. Thanks ahead!
Note that the pseudo-code in the paper is not particularly good.
The actual implementation differs in a few places. I used the C++ code to guide me.

Searching is also slow but the knn_search uses the same expensive search_layer function
that takes so much time during creation.
If we speed up creation, knn_search should improve as well.

[zgornel]

Short update:

• Read the paper and (quickly) the implementation;
• The Neighbor and NeighborSet types as well as associated functions seems pretty optimal to me, there is very little to improve there (only basic array operations); had doubts about searchsortedfirst, turns out it is very fast (takes nanoseconds even for large arrays);
• Distances.jl is pretty fast ... maybe precompute distances ?(i.e. by pairwise); will think about it for search_layer
• Will look more into the code.

Could you provide some estimate on the magnitude of the performance discrepancy between the C++ implementation and this one ?
I believe you are correct for @simd, not really sure. The answer may lie here
So far I can see one big risk: the julia implementation is slower because of many small slightly slower operations (insert/pop/delete/distance calculation etc).
Best,

[JonasIsensee]

I've also thought about the pairwise precomputation of distances.
As far as I understand though, this will only be fast if the data matrix is a proper/contiguous one.
This would then likely involve a lot of copying of data.
Could be worth it, depending on the speedup in distance computation.

Precomputing all distances in the beginning seems tempting for small datasets but this
is obviously O(N^2) so not an option.

I've looked at the video you linked and watched the section on @simd but I did not find anything helpful for us.

[JonasIsensee]

Concerning the runtime:
The values in brackets refer to the C++ implementation with num_threads=1.
With multithreading enabled I get another factor ×2.5 in speed. (i7 2-core)

efConstruction = 200
M = 16
metric = Euclidean()
ef = 50
k = 10
num_elements=10000 
dim=32
 # Float64
     creation ~3.6s  ( ~1.6s)
     search ~1.3s (~0.5s)
# Float32
    creation ~2.8s (~1.6s)
    search ~1s (~0.4s)

dim = 128
# Float64
    creation ~8.9s (~3.8s)
    search ~3.8s (~1.1s)
# Float32
    creation ~5s (~3.6s)
    search ~2.3s (~1.1s)

dim=127
# Float64 
    creation ~9.9s (~9.9s)
    search ~4.3s (~2.5s)

# Float32
    creation ~7.5s (~10.3s)
    search ~3s (~2.5s)

dim=1024
# Float32
    creation ~17s (19.5s)
    search ~7s (5.4s)

[zgornel]

Indeed, for @simd there is no word in the video on the range of the loop; moreover only multiple of 2 ranges are used :) . I cannot tell from simdloop.jl either. It must be that LLVM handles that automatically (internally, the CPU most probably adds zeros for the reminder of unused summation units and that may be a nop...). Will check if any explicit optimization (i.e. @simd, @fastmath, @inbounds) will improve performance.
Thanks for the running times updates.

[JonasIsensee]

I added another runtime for dim=127.
→ Interesting behaviour.
Both implementations take a performance hit. But the C++ version is hit harder..

[zgornel]

@JonasIsensee Hi,
Looked into using Distances.jl's colwise function to compute distances in the search_layer function. Short answer: not feasible (at all). A few arguments not to do it:

• colwise requires a matrix as input which would have to be created (from the neighbor indices)
• the neighbors have to be collected
• one could use an additional matrix for the data (along with the Vector{Vector}) however the rest of the modifications needed for the approach to work make the whole thing slower (i.e. insert points)
• even if it would be faster, the code expressiveness and similarity with the paper would be significantly decreased.

Regarding the performance macros i.e. @fastmath, @simd, there are not many places where these can be added; Most small functions are inlined as well.
All being said, its a pretty efficient implementation. Any further optimization (at least on my side) will most probably take quite a lot of time...
If you have any ideeas on improving the performance/ reducing complexity, I'd be happy to try them out.
Best,

[JonasIsensee]

Hi @zgornel ,
thanks for looking into this!

I'm wondering if one of our problems could be some kind of overhead of our distance calculations.
I added another timing for dim = 1024 where the graph creation
in the julia implementation actually got a little faster than the C++ version.
At least to me this seems like a plausible explanation.

Maybe this could be done away with by using something like SizedArray
for the vector types to make sure everything gets compiled down all the way through.
( i.e. too avoid that julia recomputes the optimal simd routine each time.
This is just guessing though. I have no idea if that actually happens every time)

[JoshChristie]

Checkout StaticArrays.jl which may help here for smaller matrix calculations.

[zgornel]

@JoshChristie worth hacking this one, will look into it.

[virtualgraham]

How do these benchmarks compare to https://github.com/KristofferC/NearestNeighbors.jl?

[zgornel]

Should be much faster (several orders of magnitude) for retrieval; indexing could be slower but for real-world usage most contexts it should be an issue.

[JonasIsensee]

How do these benchmarks compare to https://github.com/KristofferC/NearestNeighbors.jl?

This depends on the dimension of your data.
NearestNeighbors.jl is damn fast for low dimensions (< 10) but runtime increases rapidly
(exponentially !?) with increasing dimension.

This algorithm is different - it is data/metric agnostic
because it builds upon a graph structure.
This is somewhat of an overhead but it is independent of data dimension.
Therefore the runtime only depends on metric evaluations.
(which typically grows linearly for large dimensions)