The speedup is definitely off. See https://doi.org/10.1109/ICSPCS.2017.8270510 for an alternative implementation (which is part of Armadillo: https://arma.sourceforge.net/docs.html#kmeans).

arXiv version: https://arxiv.org/abs/1707.09094

Have you checked whether the distance computation itself is parallelized? If it's not, simply parallelizing that instead might give significantly more speedup. (Note that mlpack doesn't require OpenMP to be available or enabled, so going beyond #pragma omp parallel for can be a little bit tricky, if you want to go the route of the paper @conradsnicta linked to, but it can definitely be done.)

@conradsnicta @rcurtin I experimented with various OpenMP approaches. The most effective method is dividing the calculations into blocks.
On my machine, the benchmark results are:

~/gsoc/mlpack-benchmarks/src/kmeans$ ./naive_kmeans
Running benchmark for naive KMeans...
Benchmark for naive KMeans completed.
Total time: 5.10537 seconds

The problem is that blockSize needs to be adjusted for the specific system (That's why the memory tests currently fail). I could add it as another parameter, but I'm not sure if that's a valid/best practise. Is there maybe a solution where I could set that size dynamically?

@MarkFischinger I don't think hardcoding the block size is the way to go. It would be better to divide the available data into segments based on the number of available threads, which is determined at run-time.

The nominal segment size can be simply number_of_data_vectors / number_of_threads in integer form. The last segment will need to have a possibly different size than the nominal size, to take into account that the division may not be clean.

As an example, say we have 123 vectors and 4 threads, resulting in a nominal segment size of 30. The first 3 segments would each have a size of 30, thereby taking up 90 vectors in total. The last segment size would then be 123 - 90 = 33.

A bunch of corner cases would also need to be taken care of, such as ensuring a minimum number of vectors per thread, and dealing with the possibility of having more available threads than vectors.

I fixed the error and implemented a new minVectorsPerThread variable. Here are the new benchmarks:

mark@mark:~/gsoc/mlpack-benchmarks/src/kmeans$ ./naive_kmeans
Running benchmark for naive KMeans...
Benchmark for naive KMeans completed.
Total time: 7.37996 seconds

@rcurtin I updated the code following the feedback and optimized the distance calculation and OpenMP parallelization. Below are the additional benchmarks comparing the performance before and after these changes:

Before Optimization:

Running benchmark for naive KMeans...
Benchmark for naive KMeans completed.
Total time: 11.494 seconds
Running test for KMeans...
Assignments:         1        1        0        0

Centroids:
   5.5000   1.5000
   5.5000   1.5000

Test passed successfully!

After Optimization:

Running benchmark for naive KMeans...
Benchmark for naive KMeans completed.
Total time: 6.35895 seconds
Running test for KMeans...
Assignments:         1        1        0        0

Centroids:
   5.5000   1.5000
   5.5000   1.5000

Test passed successfully!

I went back to the optimized Euclidean distance formula, precomputing squared norms for both centroids and data points. This cuts down on redundant calculations, particularly when dealing with high-dimensional data. I also improved the OpenMP parallelization by making the parallel loops more fine-grained. Using a static schedule helps minimize synchronization overhead, leading to better thread utilization and faster execution. Finally, I focused on memory and cache efficiency by reusing precomputed values and that the data access is more cache friendly.

The reason I am hammering so hard (now for I believe the fifth time) about the distance calculation is because it is not always the right thing to simply observe that something is faster and then do it. The real question, when observing something like this, is understanding why the other thing is slower. (This is going to be a long comment...)

Starting from the basics, we are talking about the following two approaches:

• Sum (a[i] - b[i])^2 for each dimension i
• Precompute (a[i] * a[i]) and (b[i] * b[i]), then compute (a_norm[i] + b_norm[i] - 2 * a[i] * b[i]).

From a standpoint of what a computer should do, the first approach should be faster for several reasons:

• Only one pass over each a and b (nicer on memory)
• One subtract and multiply per dimension, vs. three multiplies and two add/sub operations
• No need for auxiliary memory to be allocated

My intuition is that memory latency will be the biggest issue here.

So the question, from my perspective, is not "which is faster", but instead "why is the first one not faster", if that's what you are observing. One key observation is that the dot-based approach will use OpenBLAS (which is very heavily optimized, and also uses OpenMP for parallelization, but, I think that should not come into play here). I would imagine that Armadillo has the same optimizations going on, but perhaps the compiler is not emitting the code we want. I therefore wrote up the following code:

#include <armadillo>

double compute_ed(const arma::mat& set1, const arma::mat& set2, const size_t points)
{
  double sum = 0.0;
  for (size_t i = 0; i < points; ++i)
  {
    for (size_t j = 0; j < points; ++j)
    {
      sum += arma::norm(set1.col(i) - set2.col(j), 2);
    }
  }
  return sum;
}

double compute_ed2(const arma::mat& set1, const arma::mat& set2, const size_t points)
{
  double sum = 0.0;
  const size_t n_rows = set1.n_rows;
  for (size_t i = 0; i < points; ++i)
  {
    for (size_t j = 0; j < points; ++j)
    {
      double point_sum = 0.0;
      for (size_t k = 0; k < set1.n_rows; k += 4)
      {
        // take 4 elements at a time to try and "trick" the vectorizer into actually working
        const double a1 = set1.mem[k + i * n_rows];
        const double a2 = set1.mem[k + 1 + i * n_rows];
        const double a3 = set1.mem[k + 2 + i * n_rows];
        const double a4 = set1.mem[k + 3 + i * n_rows];

        const double b1 = set2.mem[k + j * n_rows];
        const double b2 = set2.mem[k + 1 + j * n_rows];
        const double b3 = set2.mem[k + 2 + j * n_rows];
        const double b4 = set2.mem[k + 3 + j * n_rows];

        const double ps1 = (a1 - b1) * (a1 - b1);
        const double ps2 = (a2 - b2) * (a2 - b2);
        const double ps3 = (a3 - b3) * (a3 - b3);
        const double ps4 = (a4 - b4) * (a4 - b4);

        point_sum += ps1 + ps2 + ps3 + ps4;
      }
      sum += std::sqrt(point_sum);
    }
  }
  return sum;
}

double compute_ip(const arma::mat& set1, const arma::mat& set2, const arma::vec& set1_norms, const arma::vec& set2_norms, const size_t points)
{
  double sum = 0.0;
  for (size_t i = 0; i < points; ++i)
  {
    for (size_t j = 0; j < points; ++j)
    {
      sum += std::sqrt(set1_norms[i] + set2_norms[j] - 2 * arma::dot(set1.col(i), set2.col(j)));
    }
  }
  return sum;
}

int main()
{
  const size_t dims = 100;
  const size_t points = 1000;

  arma::mat set1(dims, points, arma::fill::randu);
  arma::mat set2(dims, points, arma::fill::randu);

  arma::wall_clock c;

  // Compute distances the EuclideanDistance way.
  c.tic();
  double sum = compute_ed(set1, set2, points);
  const double euclidean_time = c.toc();
  std::cout << "norm()-based calculation: " << euclidean_time << "s (result "
      << sum << ")." << std::endl;

  c.tic();
  sum = compute_ed2(set1, set2, points);
  const double euclidean_2_time = c.toc();
  std::cout << "norm()-based manual calculation: " << euclidean_2_time << "s (result "
      << sum << ")." << std::endl;

  // Now compute the norm-based way.
  c.tic();
  arma::vec set1_norms(points, arma::fill::none);
  arma::vec set2_norms(points, arma::fill::none);

  for (size_t i = 0; i < points; ++i)
  {
    set1_norms[i] = arma::dot(set1.col(i), set1.col(i));
  }
  for (size_t i = 0; i < points; ++i)
  {
    set2_norms[i] = arma::dot(set2.col(i), set2.col(i));
  }
  const double selfnorm_time = c.toc();
  std::cout << "dot()-based norm preparation time: " << selfnorm_time << "s."
      << std::endl;

  c.tic();
  sum = compute_ip(set1, set2, set1_norms, set2_norms, points);
  const double norm_time = c.toc();
  std::cout << "dot()-based calculation: " << norm_time << "s (result " << sum
      << ")." << std::endl;
}

Here I implemented compute_ed2() because I noticed that the Armadillo code in compute_ed() was not being vectorized properly to use SIMD instructions. When I compile with g++ -std=c++17 -march=native -DARMA_NO_DEBUG -O3 -o dot_test dot_test.cpp -larmadillo and run, I see:

$ ./dot_test 
norm()-based calculation: 0.0977767s (result 4.07441e+06).
norm()-based manual calculation: 0.0344042s (result 4.07441e+06).
dot()-based norm preparation time: 0.000114631s.
dot()-based calculation: 0.0290754s (result 4.07441e+06).

So you can see that I am getting compute_ed2() close-ish (like within 20% roughly). The big difference between compute_ed1() (which uses arma::norm() directly) and compute_ed2(), I believe, is that the vectorizer fails to properly vectorize the arma::norm() loop. Glancing into the code, I am wondering whether a temporary Armadillo object is being made for the expression set1.col(i) - set2.col(j), but I have not dug deep enough to be sure either way.

Digging deeper would involve compiling with -fopt-info-vec-missed and other similar things to see what exactly is going on. But where I am going with this is that the original distance computation strategy should be at least equivalently fast, and if it is not, we need to put in time to understand why instead of switching it out.

CC: @conradsnicta, would be interested to see what you think. I suppose we could be using norm() wrong or there is something we should use instead. I don't mind putting in some time to debug and optimize this (since distance computations are a huge part of mlpack!) but I think it would be useful to get your quick thoughts on the problem first.

@rcurtin @MarkFischinger Interim hazy perspective below. I'm currently recovering from surgery, so brain is working at reduced capacity.

For simple expressions where the memory is directly accessible in a linear column-major order (so excluding subrows and any compound expressions), the functions arma::dot() and arma::norm() generally just call the corresponding BLAS functions for common cases. This may result in using accelerated OpenBLAS implementations. For arma::norm() there is also an extra check for underflows and overflows, with a fallback to a slow but robust algorithm.

For compound expressions and subrows, arma::dot() and arma::norm() will use built-in algorithms to avoid making copies of data. In some cases these algorithms can be vectorised by the compiler, though the speed difference is around 20% in my (limited) observations on an ancient laptop. The speedup may be higher with CPUs that handle AVX-512.

There is a bunch of caveats here:

• Various compilers will differ in what they're willing to vectorise, in order to stick with the IEEE floating point standard. For GCC, more vectorisation can occur if the -ffinite-math-only flag is used, but this can come with the possible side-effect of not treating NaN and Inf values properly. Haven't fully checked how picky clang is.

• AVX instructions generally don't get used with GCC (and possibly clang) unless the -march=native flag is used.

• For GCC, need to use the -O3 optimisation level for vectorisation to be fully enabled. I think clang enables vectorisation at the -O2 level, though its vectoriser may not be as advanced as GCC's vectoriser at -O3.

Other observations:

• For the case of the k-means algorithm, using arma::norm() for the distance function is doable, but there is an unnecessary call to std::sqrt() which will cost a bit of time.

• I suggest comparing the speed of the modified k-means implementation against arma::kmeans(). This function is automatically multi-threaded if OpenMP is enabled via the -fopenmp option.
  https://arma.sourceforge.net/docs.html#kmeans

@conradsnicta and @rcurtin, thanks for your feedback. I’ve been digging deeper into the performance of the sum of squared differences versus the dot-product based approach. Despite initial expectations that the sum of squared differences would be faster, my benchmarks show that both methods perform almost identically, prompting further investigation.

Here are the initial benchmark results:

Sum of squared differences time: 1.02083 seconds, result: 1.65892e+08
Dot-product based time: 1.037 seconds, result: 1.65892e+08

After applying the -fopt-info-vec flag, both methods were optimized with SIMD instructions using 16-byte vectors from the AVX instruction set. The updated benchmarks were:

Sum of squared differences time: 1.04866 seconds, result: 1.66398e+08
Dot-product based time: 1.03961 seconds, result: 1.66398e+08

Enabling the -ffast-math flag led to a significant reduction in execution times:

Sum of squared differences time: 0.272451 seconds, result: 1.67351e+08
Dot-product based time: 0.253304 seconds, result: 1.67351e+08

These results show how aggressive compiler optimizations impact both methods, bringing their performance closer together.

I am wondering whether a temporary Armadillo object is being made for the expression set1.col(i) - set2.col(j)

I also considered whether temporary Armadillo objects during operations like set1.col(i) - set2.col(j) could be affecting performance. However, due to the use of Proxy<T1> P(X);, which efficiently manages references, this was ruled out as a factor.

For further comparison, I benchmarked our implementation against arma::kmeans().

Benchmark for mlpack KMeans completed.
Total time: 1.09256 seconds
Benchmark for Armadillo KMeans completed.
Total time: 0.00981316 seconds

I'm a bit unsure why the performance diffenerence is so huge. Armadillo is more optimized, but I'm looking into it.

When reviewing the assembly code, I observed that both the sum_of_squared_differences and dot_product_distance functions are using SIMD instructions, particularly with AVX (Advanced Vector Extensions) and FMA operations like vfmadd231pd.
However, the performance gap between the mlpack implementation and arma::kmeans() suggests there’s considerable room for improvement.

Currently, I’m focusing on improving the arma::norm function to see if we can apply similar optimizations.

It's not possible to directly compare mlpack's and Armadillo's k-means implementations without further digging. They will run for different numbers of iterations, to different tolerance measures, etc., etc... the way to compare them with each other would be to time how long a single iteration takes, or how long a specific number of iterations takes. (That's not to say that the iteration or termination strategies in mlpack couldn't be improved. But here we are considering only the distance computations, so ensuring they are equal across both implementations is needed.)

@MarkFischinger, looking at this PR, it seems to me that it is getting further from the objective. For the sake of GSoC, please could you do the following:

1. Revert the modifications related to the distance calculation
2. keep the original update regarding the OpenMP
3. Do not compare mlpack and armadillo implementations since they are not relevant.

Please do all of the above requests by the end of today, if you are interested in optimising the distance calculation let us do that in a different PR and merge this one with OpenMP modification only

Thank you very much.

I do think that we need to get to the bottom of why arma::norm() isn't giving good performance in our case, but for the sake of getting something merged, I think just OpenMP-izing the distance calculations (and leaving the optimization of the individual distance calculations for later) is sufficient. We can return to the arma::norm() issue elsewhere---and working that out will have significant benefits throughout the library because of how widely distance calculations are used. 👍

There are a lot of changes here that aren't actually functional changes---do you mind going through and reducing the changes to only the actual functionality changes? Thanks 👍

@rcurtin @shrit I added all the important comments back in, replaced .unsafe_col() with .col() like you suggested, fixed the line issue at the end and removed the arma::zeros, because of this comment.

Running the benchmark again using this benchmark script and this visualization script, I generated this graphic:

Thanks for the hard work @MarkFischinger! I opened #3789 to follow up on the EuclideanDistance efficiency issue. We can continue that there.

I think I need to fix the documentation build so that it issues warnings and not errors---in this case, looks like the CS department's website for Boston College has gone offline, causing that build to fail. Oh well.