I'm working on a small fun project on the side for a numerical method I've been working a bit with. Roughly, the computational problem I have to solve is the following: Assume you have a collection of $N$ points, call this collection $X$. For each $x_i \in X$, the following steps occur:

1: Find $k$ nearest neighbors to $x_i$, call this set $X_i$.

2: Build a $k \times k$ matrix out of these points with some simple arithmetic. Call this matrix $A^i$. $A^i(j,k)$ depends only on points $y_j,y_k \in X_i$.

3: Build a vector $b^i$ who is 1 only in the row corresponding to the originally chosen point $x^i$.

4: Solve the linear system $A^if^i =b^i$ with whatever method is best.

5: Return the $k \times 1$ vector $f^i$ from Step 4 as well as the indices which map points in $X^i$ back to $X$ (since $X^i$ is a subcollection of $X$).

So we have those 5 steps to run through for each point $x_i \in X$. So, if $N = 30,000$, for example, then we have to run through Steps 1-5 30,000 times. Each iteration is independent: there needs to be no communication between running iteration $i$ and iteration $j$. They all rely on the same common dataset, which is $X$.

I've done this in Matlab before and it was pretty easy to code up, but even with vectorizing as much as I could, it got a bit slow when N = 30,000 and k = 400. I decided as a good lesson for myself, to try to write it in C++.

My thoughts on the C++ code were the following:

1: Step 1 requires computing nearest neighbors. Perhaps it would be helpful to use a tree structure to keep in mind what points are nearby. I used Boost's spatial library to set up an R Tree and I query that for step 1 to build the set $X_i$. Please see http://www.boost.org/doc/libs/1_58_0/libs/geometry/doc/html/geometry/spatial_indexes/introduction.html

2: Since i needed linear algebra, but my matrices were "small", I figured something like Trilinos or Petsc would be overkill, so I went with Armadillo. I use armadillo to build the matrices $A^i$, the vectors, and handle the solve step.

3: I'm not very experience with OOP, so I may have made a major faux pas here. I made a class Assembler that has the dataset $X$ and the Rtree as member variables (it builds them), and then it has a method I can call that will instantiate an object L, which L then runs steps 1-5 and stores the solution $f^i$ and the indices as member variables.

4: I have what looks to me to be an embarassingly parallel problem, so I call my Assembler class's method to build up the individual objects L and stick that in a for loop for i=1 to i= N. I then just stuck an OpenMP pragma around the for loop. This gave me pretty good speedups (better than Matlab by a lot, depending on what machine I ran it on).

My questions are the following: Am I on the right track? Anyone have any suggestions or advice on how to improve things? Any suggestions on speed-ups or improved parallelism? If it would be helpful, I can get a link to the code so people can take a look at it.


EDIT: Another question I forgot to ask: each iteration of building object L_i runts through Steps 1-5, but at the end of the day, all I want is to store (for later use) the coefficient vector $f^i$ and the index vector mapping $X_i$ back to $X$ (both of which are member variables object $L_i$ finds during Steps 1-5). Naively, I'm inclined to write these to a csv file or binary file or something, and have one file per object $L_i$, but I'm not exactly sure if that's smart. In MATLAB, I had a sparse matrix $S$ where column $i$ of $S$ was the vector $f^i$ indexed properly with the map that sends $X_i \to X$. Since $k << N$, $S$ was a sparse matrix. This got super slow in MATLAB, as I guess I wasn't allocating $S$ properly and it was extending itself or something.

  • $\begingroup$ A minor point: when solving $Ax=b$ with LU factorization and $b$ has zeros at the beginning, the initial steps of forward substitution $L^{-1}b$ are redundant and can be skipped, which saves a bit of time. $\endgroup$
    – Kirill
    Dec 8, 2015 at 18:44

2 Answers 2


The linear algebra part of the problem looks like something that would be ideal for use on GPUs; every matrix has the same size and there are a whole bunch of them. You can use the routines available in the ViennaCL for this. ViennaCL has backends for OpenMP, OpenCL and CUDA, so you can, in principle, use it on a multicore CPU or on any GPU; you're not restricted to just NVIDIA hardware as you would be with a CUDA library like cuBLAS. There's also clBLAS and any number of alternatives.

If you decide to use GPUs, you might also be able to build the matrices on the GPU directly, which will save the overhead of moving the memory from the CPU to the GPU. Since you have 30000 matrices of size 400 x 400, this would be a huge savings.

If you decide not to go with GPUs, you could also try using a library like ATLAS together with Armadillo, which may provide faster dense linear algebra kernels than the default ones that Armadillo uses. There's more information about using optimized linear algebra libraries on the Armadillo FAQ.

You also mentioned using a class to assemble the matrices for the problem and were concerned with whether or not this was a good idea. I doubt this is going to be much of a concern. The performance penalties you run into with object-oriented programming are usually a result of dynamic dispatch, which requires a level of pointer indirection. However, this will be minuscule compared to the floating point operations in building a 400 x 400 matrix. You may want to investigate other spatial data structures than R-trees, like using a Z-ordering and throwing out the tree data structure altogether. My suspicion is that this will, again, pale in comparison to the linear algebra operations.

This brings me to my final point: a good profiling tool (scalasca, TAU, ompP) will tell you far more about where your code spends the most time, and therefore how to optimize it, than some guy on the internet will.

  • $\begingroup$ Thanks for the response! So, right now I've got Armadillo linked against OpenBLAS on one of my machines, and on another machine it's using ACML for BLAS and LAPACK;that definitely gave a very nice speedup. Using GPUs was my first instinct as each matrix can be built up, on GPU, using the same dataset $X$. Unfortunately, I didn't find much NVidia support for solving multiple matrices at a time on GPU. CUDA 7 has a batch QR for multiple sparse matrices, but no such luck for many systems of matrices already on GPU (or I'm misunderstanding their docs). I will definitely start profiling! $\endgroup$
    – user35959
    Dec 6, 2015 at 16:44
  • 1
    $\begingroup$ You may have to write some small amount of CUDA/OpenCL code to solve a whole batch of systems; the Nvidia libraries will provide you with just the routines to solving a single dense system. AMD has deprecated ACML in favor of open-source libraries; these might have improvements over ACML, so you could see if Armadillo can use them too. $\endgroup$ Dec 6, 2015 at 22:57
  • For finding the nearest neighbors quickly, a KD-Tree would probably be faster than a R-Tree: a R-Tree supports dynamic updates, that you do not need (and you pay a price for it). Some implementations of KD-Tree are avaible in David Mount's ANN library [1] and in my Geogram library [2]. Note that ANN is not thread-safe (but it is quite easy to make it tread-safe by changing a little bit its source-code / replacing some global variables by function call parameters). Geogram is thread-safe. There is also a FLANN library [3] that is thread-safe (but with this one I did not manage to have exactly the nearest neighbors, there was always an approximation, but maybe I did not use it correctly).

  • Another word on nearest-neighbors: if you do not need exactly the k nearest neighbors but can do some approximations, then you may gain a little bit, see the documentation of ANN and FLANN.

  • With a size of 400x400, maybe using a standard dense linear solver will be the fastest, such as those in LAPACK, especially if you find a version of LAPACK optimized for your processor, using vector instructions (AVX2). GPU may be faster, but it can be a pain to setup and to program. If you want something more modern than LAPACK, you may try also eigen [4], it has many different solvers.

To summarize, to answer your question, I would do a big #pragma omp parallel for loop, that retreives the neighbors using a Kd-Tree and solves the linear system using dense linear algebra with AVX instructions for vectorization. In the loop, avoid doing any dynamic memory allocation, since it breaks multithreading (on most OS/c++ runtimes, there is a global lock on dynamic memory allocation/deallocation). Use instead local variables (allocated on the stack, each thread has an independent one).

Last but not least : I agree with D. Shapero's remark, before starting to optimize anything, it is very important to run a profiler: each time I did that, I was surprised, the bottleneck was absolutely not where I expected it !!

[1] https://www.cs.umd.edu/~mount/ANN/

[2] http://alice.loria.fr/software/geogram/doc/html/index.html

Geogram Kd-Tree class: http://alice.loria.fr/software/geogram/doc/html/kd__tree_8h.html

[3] http://www.cs.ubc.ca/research/flann/

[4] http://eigen.tuxfamily.org

  • $\begingroup$ Thanks for the info! I actually wanted a Kd Tree from the start, but a quick google didn't help much (I was looking around for one in Boost). Your Geogram looks great, I'll try to grab that and get it running. I also appreciate the suggestion of avoiding dynamic allocation. That hadn't occurred to me, and I do use dynamically allocated structures (std::vecetor<T>'s throughout). I'll also start profiling per your and Daniel's remark and probably post an update later. $\endgroup$
    – user35959
    Dec 6, 2015 at 16:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.