Advice on handling many "small" matrices in parallel

Question

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.

Thanks!

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.

Answer 1

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.

Answer 2

• 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